Matlab Numerical Solver: Solving Second Order differential Equation

Question

I am truly sorry that I could not provide details for the exact equation that I am working with. It is a very complicated second-order differential equation in the form similar to this:

where function a(z) ~ e(z) and g(z) are given. p is constant.

I also have the boundary conditions.

Is it possible to solve f(z) with the help of matlab?

Any suggestion would be greatly appreciated. Thanks a lot.

EDIT

Here is how I defined my function. a1 ~ g1 and fna ~ fng are defined and stored in Gdata.mat:

function xp=myfunc(t,p)
% x =  [d2f df f df2 f2 d2f2]
% xp = [df d2f f2 df2 d3f d2f2]
load GData
xp=zeros(6,1); % [f df d2f f2 df2]
               % f
fprintf('%d\n',length(xp));
fprintf('%d\n',length(p));
xp(1) = x(2);  % df
xp(2) = x(1);  % d2f
               % f2
xp(4) = x(4);  % df2
xp(6) = x(6);  % d2f2
xp(5) = (...
        b1(t)*p(3) + b(t)*p(2) + ...
        c1(t)*p(3)^3 + 3*fnc(t)*p(3)^2*p(2) + ...
        d1(t)*p(3)^5 + 5*fnd(t)*p(3)^4*p(2) + ...
        e1(t)*p(3)^7 + 7*fne(t)*x(3)^6*p(2) - ...
        f1(t)*p(2)*p(3) + f1(t)*p(1)*p(3) + f1(t)*p(2)^2 - ...
        g1(t)*p(4) - fng(t)*p(6) + ...
        q*p(2) - a1(t)*p(1)...
        ) * 1/(fna(t));

then I called:

[TEMP,POL] = ode45('odesolver',[0,1],[0,0,0,0,0,0]);

EDIT2

function dp=odesolver(t,p)
% dp = [df d2f d3f]
syms x;
load BData;
load GData;
dp=zeros(3,1); % [f df d2f]

A = interp1(t_data,At,t);
B = interp1(t_data,Bt,t);
C = interp1(t_data,Ct,t);
D = interp1(t_data,Dt,t);
E = interp1(t_data,Et,t);
F = interp1(t_data,Ft,t);
G = interp1(t_data,Gt,t);

A1 = interp1(t_data,A1t,t);
B1 = interp1(t_data,B1t,t);
C1 = interp1(t_data,C1t,t);
D1 = interp1(t_data,D1t,t);
E1 = interp1(t_data,E1t,t);
F1 = interp1(t_data,F1t,t);
G1 = interp1(t_data,G1t,t);

dp(1) = p(2); % f'
dp(2) = p(3); % f''
dp(3) = (...
        B1*p(3) + B*p(2) + ...
        C1*p(3)^3 + 3*C*p(3)^2*p(2) + ...
        D1*p(3)^5 + 5*D*p(3)^4*p(2) + ...
        E1*p(3)^7 + 7*E*p(3)^6*p(2) - ...
        F1*p(2)*p(3) + F*p(1)*p(3) + F*p(2)^2 - ...
        2*G1*p(2)*p(1) + 2*G*p(3)*p(1) + 2*G*p(2)*p(2) + ...
        q*p(2) - A1*p(1)) * 1/(A);

Answer 1

Answer based on the discussion and edited question:

There are several obstacles in using ode45 to solve your differential equation, but none of them are a showstopper:

• Second order ODE: convert into 2 first-order odes you can solver with ode45, as in this question.
• Integral term: differentiate your equation to get rid of it. You will end up with a third-order differential equation, which you need to convert into 3 first-order equations using the same technique as above
• Derivative of the square of the function to integrate: add additional states to your ode function to allow solving with ode45
• Time-dependent variables: you can't just directly use your time variable to index into your time-dependent variables. You need to do a linear interpolation to for each value of t using your time-dependent data to get the corresponding values of your variables, as in this other question.

So your code should look something like this (I have assumed that you have a common time vector for all your time-dependent data in your GData file, this may not be the case, you may have a different time vector for each variable - adjust the code if necessary; either way make sure your time-dependent data is defined on the time interval on which your are calling ode45):

function dx=myfunc(t,x)
% x =  [d2f df f df2 f2 d2f2]
% dx = [d3f d2f df d2f2 df2 d3f2]
load GData % I assume this contains b1, b, c1, etc... *and* the corresponding time vector, say t_data
dx=zeros(6,1); % [d3f d2f df d2f2 df2 d3f2]

fprintf('%d\n',length(dx));
fprintf('%d\n',length(x));
dx(1) = x(2);  % double-check!!
dx(2) = x(1);  % double-check!!
               % double-check!!
dx(4) = x(4);  % double-check!!
dx(6) = x(6);  % double-check!!

B1 = interp1(t_data,b1,t);
B = interp1(t_data,b,t);
C1 = interp1(t_data,c1,t);
% etc... for the other variables

dx(5) = (...
        B1*p(3) + B*p(2) + ...
        C1*p(3)^3 + 3*FNC*p(3)^2*p(2) + ...
        D1*p(3)^5 + 5*FND*p(3)^4*p(2) + ...
        E1*p(3)^7 + 7*FNE*x(3)^6*p(2) - ...
        F1*p(2)*p(3) + F1*p(1)*p(3) + F1*p(2)^2 - ...
        G1*p(4) - FNG*p(6) + ...
        q*p(2) - A1*p(1)...
        ) * 1/(FNA); % double-check!!