Not getting six coupled ODEs solution using bvp4c in MATLAB, Please help me

13 views (last 30 days)
I am solving six coupled first order ODEs using bvp4c in MATLAB, but it is not working...here is my code
function finalprogramme
thetaa = 0; thetab = pi/2 ;
solinit = bvpinit(linspace(thetaa,thetab,1000),[0.015 0 0 1 0 0]);
sol = bvp4c(@sixodes,@bc,solinit);
thetaint = linspace(thetaa,thetab,1000);
Sthetaint = deval(sol,thetaint);
plot(thetaint,Sthetaint([1 3],:));
function res = bc(ya,yb)
L = 70; u_star = 769.23e2; F=0.5*L;
res = [ ya(2); ya(3); ya(5); yb(2); yb(3); yb(5)-(-F/(2*L*u_star*H))];
function dydtheta = sixodes(theta,y)
R = 10; neu_star=0.3; alpha=(neu_star)/(1-2*neu_star); K=(12*(R^2)*(1-neu_star))/(H^2);
dydtheta = [y(2); ((5*K*alpha+4*alpha+3*K+3)/(1+K*(alpha+1)))*y(4)+((2*(1+2*alpha)*(K+1))/(1+K*(alpha+1)))*y(1)-R*((1-K*alpha-K)/(1+K*(alpha+1)))*y(6); y(4); ((-1)/(2*(2*alpha+1)))*((3*alpha+2)*y(2)-alpha*y(3)-alpha*R*y(5)); y(6); ((K+2)/(2))*(((-1)/(Ri))*y(2)+((1)/(R))*y(3)+y(5))];

Accepted Answer

Stephan
Stephan on 10 Jun 2020
Edited: Stephan on 10 Jun 2020
I used nested functions - for 2 of the constants you did not provide values:
% Call the outer function
finalprogramme
%% Outer function
function finalprogramme
% You missed to tell us values for these both constants - used 1 for them
H = 1;
Ri = 1;
% All other Constants
thetaa = 0;
thetab = pi/2;
R = 10;
neu_star=0.3;
L = 70;
u_star = 769.23e2;
F=0.5*L;
alpha=(neu_star)/(1-2*neu_star);
K=(12*(R^2)*(1-neu_star))/(H^2);
% Solve the problem
solinit = bvpinit(linspace(thetaa,thetab,1000),[0.015 0 0 1 0 0]);
sol = bvp4c(@sixodes,@bc,solinit);
thetaint = linspace(thetaa,thetab,1000);
Sthetaint = deval(sol,thetaint);
plot(thetaint,Sthetaint([1 3],:));
%% Inner functions
function res = bc(ya,yb)
res = [ ya(2); ya(3); ya(5); yb(2); yb(3); yb(5)-(-F/(2*L*u_star*H))];
end
function dydtheta = sixodes(~,y)
dydtheta = [y(2); ((5*K*alpha+4*alpha+3*K+3)/(1+K*(alpha+1)))*y(4)+((2*(1+2*alpha)*(K+1))/(1+K*(alpha+1)))*y(1)-R*((1-K*alpha-K)/(1+K*(alpha+1)))*y(6); y(4); ((-1)/(2*(2*alpha+1)))*((3*alpha+2)*y(2)-alpha*y(3)-alpha*R*y(5)); y(6); ((K+2)/(2))*(((-1)/(Ri))*y(2)+((1)/(R))*y(3)+y(5))];
end
end
Result (using my constants for H & Ri) is:

More Answers (0)

Categories

Find more on Function Creation in Help Center and File Exchange

Tags

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!