How do I use 'deval' or any equivalent kind in the code for shooting method using 'ode15s'

Hello,

I need to find out the value of y(3),y(5), y(7) at x=0,

There is an option 'deval' eg. deval(sol,0), which I used in bvp4c solver, but I have no idea how to use it in this particular solver i.e ode15s

Following is my code:

function shooting_method
clc;clf;clear;
global lambda gama Pr Rd Lew Nb Nt Mn
gama=1;  
Mn=1;
Rd=0.1;
Pr=10;
Nb=0.2;
Lew=10;
Nt=0.2;
lambda=0.5;
x=[1 1 1];
options= optimset('Display','iter');
x1= fsolve(@solver,x);
end
function F=solver(x)
options= odeset('RelTol',1e-8,'AbsTol',[1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8]);
[t,u]= ode15s(@equation,[0 4],[0 1 x(1) 1 x(2) 1 x(3)],options)
s=length(t);
F= [u(s,2),u(s,4),u(s,6)];
%deval(0,u)
plot(t,u(:,2),t,u(:,4),t,u(:,6));
end
function dy=equation(t,y)
global lambda gama Pr Rd Lew Nb Nt Mn
dy= zeros(7,1);
dy(1)= y(2);
dy(2)= y(3)*(y(3)^2+gama^2)/(y(3)^2+lambda*gama^2);
dy(3)= y(2)^2/3-(2*y(1)*y(3)*(y(3)^2+gama^2))/(3*(y(3)^2+lambda*gama^2))+Mn*y(2);
dy(4)= y(5);
dy(5)= -(2*Pr*y(1)*y(5))/(3*(1+Rd)) - (Nb*y(5)*y(7))/(1+Rd) - (Nt*y(5)^2)/(1+Rd);
dy(6)= y(7);
dy(7)= -((2*Lew*y(1)*y(7))/3)+ (Nt/Nb)*((2*Pr*y(1)*y(5))/(3*(1+Rd)) + (Nb*y(5)*y(7))/(1+Rd) + (Nt*y(5)^2)/(1+Rd));
end

Answer 1

[Harshavardhan]

Hi @naygarp

I see that you are trying to find the value of the variables y(3),y(4),y(5) at x=0 using ode15s as the solver. To do so, you can use the function “interp1” as shown below.

% Interpolate at x = 0  
y_at_x0 = interp1(t, u, 0, 'pchip'); 

% Extract required values
y3_x0 = y_at_x0(3);  
y5_x0 = y_at_x0(5); 
y7_x0 = y_at_x0(7);  

fprintf('y3(0) = %.5f, y5(0) = %.5f, y7(0) = %.5f\n', y3_x0, y5_x0, y7_x0);  

For more information on “interp1” you can refer to the following link :

https://www.mathworks.com/help/matlab/ref/double.interp1.html

or type the following command in a MALTAB command window for version specific documentation.

doc interp1