MATLAB: Parameter variation in BVP shooting method

bvpnumericsode45shooting method

Hi everybody!
I wrote a code that solves the following problem:
This is a 2nd order ode and the dependant variable here is rho (tilde) and the independant variable is x (tilde). K (tilde) and B (tilde) are both parameters. The boundary conditions are:
The code is:
clear all
clc
%Shooting Method for Non-Linear ODE
k=3.25e8; %Debye length - needs to be calculated.
D=4e-9; %separation distance
y0=1; %i.c
ICguess_on_target=fzero(@(x) bar_res(x), 0.25);
[x,y]=ode45(@bar_density,[0,k.*D],[y0 ICguess_on_target]);
plot(x,y(:,1));
xlabel('\kappaD');
ylabel('$\tilde{\rho}$','Interpreter','latex')
title('$\tilde{\kappa}=\tilde{B}=1$','Interpreter','latex'); 
function dydx=bar_density(x,y)
kdim=1;Bdim=1;
dydx=[y(2);2.*Bdim./kdim.*y(1)+kdim.^(-1).*log(y(1))-Bdim];
end
function r=bar_res(IC_guess)
y0=1;
yf=1;
k=3.25e8;
D=4e-9;
[x,y]=ode45(@bar_density,[0,k.*D],[y0 IC_guess]);
r=y(end,1)-yf
end
This works just fine and I'm getting great results and plots of rho(tilde) as a function of x(tilde). Now I want to run the code in a loop and vary the parameters K (tilde) and B (tilde). For instance I would like to get 6 plots of rho(tilde) as a function of x(tilde) for 6 different sets of K (tilde) and B (tilde). I have been trying to insert a for loop inside the code but I'm obviously doing it wrong.
I would appreciate some help guys.
Thanks,
Roi

Best Answer

The following works:
% Shooting Method for Non-Linear ODE
k=3.25e8; %Debye length - needs to be calculated.
D=4e-9; %separation distance
y0=1; %i.c
kdimarray = 1:0.2:2; Bdimarray = 1:0.5:3.5;
ICguess_on_target = zeros(size(kdimarray));
for i = 1:6
    kdim = kdimarray(i); Bdim = Bdimarray(i);
    ICguess_on_target(i) = fzero(@(x) bar_res(x,kdim,Bdim), 0.25);
    [x,y]=ode45(@bar_density,[0,k.*D],[y0 ICguess_on_target(i)],[], kdim, Bdim);
    plot(x,y(:,1));
    hold on
end
xlabel('\kappaD');
ylabel('$\tilde{\rho}$','Interpreter','latex')
title('$\tilde{\kappa}=\tilde{B}=1$','Interpreter','latex'); 
function dydx=bar_density(~, y, kdim, Bdim)
dydx=[y(2);2.*Bdim./kdim.*y(1)+kdim.^(-1).*log(y(1))-Bdim];
end
function r=bar_res(IC_guess,kdim,Bdim)
y0=1;
yf=1;
k=3.25e8;
D=4e-9;
[~,y]=ode45(@bar_density,[0,k.*D],[y0 IC_guess], [], kdim, Bdim);
r=y(end,1)-yf;
end
Related Question