MATLAB: How to solve a set of 4 first-order non-linear coupled ODEs

nanode

I am trying to solve a 4th order differential equation using shooting method by disintegrating the ODE into four coupled first order ODEs. I only have initial conditions. So I am trying to perturb y(1) and y(4) and use them as a parameter to achieve a target value for y(3). But all I get is NaNs on my output matrices for y(1), y(2), y(3) and y(4). Any inputs will be appreciated.
Here is my code:
 %%Constants
 A = 1e-20;          % Hamaker Constant
 R_FC = 24.597;      % Gas Constant of FC-72 (J/kg-K)
 c = 1;              % Accomodation Coefficient
 sigma = 0.01;       % Surface Tension Coefficient (N/m)
 rho_l = 1593.84;    % Liquid Density (kg/m3)
 k_l = 0.0544;       % Liquid Conductivity (W/m-K)
 mu_l = 0.0004683;   % Liquid Viscosity (kg/m2-s)
 rho_v = 11.61;      % Vapor Density (kg/m3)
 h_lv = 88000;       % Enthalpy of Vaporization (J/kg)
 Twall = 334;    % Wall Temperature (K), 10 degrees superheat
 Tsat = 329;     % Saturation Temperature of FC-72 at 1 atm (K)
 R_int = ((2 - c)/(2*c))*(Tsat^(1.5))*((2*pi*R_FC)^(0.5))/(rho_v*(h_lv^2));   
 delta_ad = (A/(rho_l*h_lv*((Twall/Tsat) - 1)))^(1/3);
 e1 = delta_ad/1000;
 deltaP_ad = A/((delta_ad)^3);
 deltaP_ad_corr = A/((delta_ad+e1)^3);
 eQ = 1;
 %%Function
 tpcl = @(xi,y) [y(2); (1/sigma)*((1+(y(2)*y(2)))^(1.5))*(y(3)-(A/(y(1)^3))); (3*mu_l/(rho_l*h_lv))*(-y(4)/(y(1)^3)); (Twall-Tsat-(Tsat*y(3)/(rho_l*h_lv)))/((y(1)/k_l)+R_int)];
xi_end = 5e-7;
N = 1000;
[xi,y1] = rk4(tpcl,[0 xi_end],[delta_ad 0 deltaP_ad 0],N);
[xi,y2] = rk4(tpcl,[0 xi_end],[delta_ad+e1 0 deltaP_ad_corr eQ],N);
delta_a = delta_ad;
delta_b = delta_ad + e1;
Q_a = 0;
Q_b = eQ;
slope_a = y1(2,end);
slope_b = y2(2,end);
 %%Shooting Method
 tol = 1e-5;
 iter = 10;
 slope_target = tan(pi/18);
 for i = 1:iter
    Q_new = Q_a + (Q_b-Q_a)/(slope_b-slope_a)*(slope_target-slope_a);
    delta_new = delta_a + (delta_b-delta_a)/(slope_b-slope_a)*(slope_target-slope_a);
    deltaP_new = A/((delta_new)^3);
    [xi,y] = rk4(tpcl,[0 xi_end],[delta_new 0 deltaP_new 0],N);
    fprintf('iter:%2d, delta(0)=%17.15f, slope(xi_end)=%17.15f\n',i,delta_new,y(2,end));
    if (abs(y(2,end)-slope_target) <= tol)
        break;
    end
    Q_a = Q_b;
    Q_b = Q_new;
    delta_a = delta_b;
    delta_b = delta_new;
    slope_a = slope_b;
    slope_b = y(2,end);
end

Best Answer

Anonymous functions are hard to debug. Better define the ODE as a function and use http://www.mathworks.de/matlabcentral/answers/1971 to define the constants on demand.
Then you can use the debugger to step through the program line by line until you find the source of the NaN's by your own. This is more efficient to let all forum users debug the code for you, especially when the function rk4() is not known.
Related Question