MATLAB: Bvp4c singular jacobian error

bvp4csingular jacobian error

Hello,
I am trying to solve a boundary value problem. But I have faced with " Singular Jacobian Error". Does anybody know how I can solve it?
It doesn't work for Nmf>4 !!!
clc; close all; clear all;
format long e
solinit = bvpinit(x,@xinitfun);
sol= bvp4c(@odefun,@bcfun,solinit);
x= linspace(-1,1,100)  
function xinit= xinitfun(x)
xinit=rand((2*Nmf+1)*6 * 2,1);
function dydx = odefun(x , y)  
global Nmf Re alpha Ra Pr NR
dydx = zeros((2 * Nmf + 1) * 2 * NR,1);
dydx1= zeros((2 * Nmf + 1) * NR,1);
for ii = 1 : 2 * Nmf + 1
      RI = (ii - 1) * NR;
      II = (ii - 1) * NR + (2 * Nmf +1) * NR;
      n = ii - Nmf - 1;
      na = n * alpha;
      teta_f2 = 0;
      NLT1 = 0;
      NLT2 = 0;
      for jj = 1 : 2 * Nmf + 1
          m = jj - Nmf - 1;
          ma = m * alpha;
          nma =(n - m) * alpha;
          teta_f = 0;
          diff_teta_f = 0;
          if abs(n - m) <= Nmf
              if abs(n - m) == 1
                  teta_f = (-1 / 8 * Pr) * sinh(nma .* x) ./ sinh(nma) + (1 / 8 * Pr) * cosh(nma .* x) ./ cosh(nma);
                  diff_teta_f = (nma * sinh(nma .* x)) ./ (8 * Pr * cosh(nma)) - (nma * cosh(nma .* x)) ./ (8 * Pr * sinh(nma));
              elseif (n - m) == 0
                  teta_f = 0.5 .* (1 - x) / Pr;
                  diff_teta_f = 0.5 / Pr;
              end
              NLT1 = NLT1 + (y(2 + RI) + y(2 + II) * i) * ma .* ((y(3 + RI) + y(3 + II) * i) - ma ^ 2 .* (y(1 + RI) + y(1 + II) * i))...
                  - nma .* (y(1 + RI) + y(1 + II) * i) .* ((y(4 + RI) + y(4 + II) * i) - ma ^ 2 .* (y(2 + RI) + y(2 + II) * i));
              NLT2 = NLT2 + nma .* (y(2 + RI) + y(2 + II) * i) .* (teta_f + Pr .* (y(5 + RI) + y(5 + II) * i)) - ma .*...
                  (y(1 + RI) + y(1 + II) * i) .* (diff_teta_f + Pr .*( y(6 + RI) + y(6 + II) * i));
          end
      end
      if abs(n) == 1
          teta_f2 = (-1 / 8 * Pr) * sinh(na .* x) ./ sinh(na) + (1 / 8 * Pr) * cosh(na .* x) ./ cosh(na);
      elseif n == 0
          teta_f2 = 0.5 .* (1 - x)/ Pr;
      end
      dydx1((ii - 1) * NR + 1 : ii * NR,1) = [y(2 + RI) + y(2 + II) * i
          y(3 + RI) + y(3 + II) * i
          y(4 + RI) + y(4 + II) * i
          2 * (na) ^ 2 .* (y(3 + RI) + y(3 + II) * i) - na ^ 4 .* (y(1 + RI) + y(1 + II) * i) + (na * Re .* (1 - x .^ 2) .* ...
          (y(3 + RI) + y(3 + II) * i - na ^ 2 .* (y(1 + RI) + y(1 + II) * i)) - na * Re * (-2) .* (y(1 + RI) + y(1 + II) * i) + na * Ra .* (y(5 + RI) + y(5 + II) * i) + (Ra/Pr) * na .* teta_f2 + NLT1) * i
          y(6 + RI) + y(6 + II) * i
          na ^2 .* (y(5 + RI) + y(5 + II) * i) + (na * Re * (1 - x .^ 2) .* (Pr .* (y(5 + RI) + y(5 + II) * i) + teta_f2 ) + NLT2) * i];
end
dydx=[real(dydx1);imag(dydx1)];
    function res=bcfun(ya,yb)
    global Nmf NR
res=[ya(1 : NR :  (2 * Nmf + 1) * 2 * NR)
     ya(2 : NR :  (2 * Nmf + 1) * 2 * NR)
     ya(5 : NR :  (2 * Nmf + 1) * 2 * NR)
     yb(1 : NR :  (2 * Nmf + 1) * 2 * NR)
     yb(2 : NR :  (2 * Nmf + 1) * 2 * NR)
     yb(5 : NR :  (2 * Nmf + 1) * 2 * NR)];

Best Answer

Related Question