# Error in ODE45, must return a column vector

54 views (last 30 days)
Brad Scott on 23 Feb 2021
Commented: Adam Danz on 26 Feb 2021
%%Can find the cbar solution using quadratic formula
%%Constants
D = (10^(-4))/(1-10^(-4));
gamma = 0.4;
M = 0.2;
z = linspace(0,1,1000);
mcQ = 10^5; %mathcal Q
n=2;
w=0.5;
a = 1;
b = D+gamma*z-M;
c = -D*M;
%%Quadratic formula, c is negative, so for real solution we just use + root
cbar = 0.5*(-b+sqrt(b.^2-4*a*c));
%Q from definition
Qbar = gamma*z + cbar - M;
%%Test solutions, avg error of -2.67*10^(-19) approx 0
test_c_Q=mean(cbar.*(D+Qbar) - D*M);
%%Find phi using root finder
for Qs=Qbar
root_est = (Qs/mcQ)^(1/n);
p(Qs==Qbar) = fzero(@(p) fphi(p,n,mcQ,Qs), root_est);
end
%%Test phi value, avg error of -6.46*10^(-17) approx 0
test_p=mean(p+mcQ*p.^(n).*(1-p).^2 - Qbar);
%%Setup to solve ODE
dQdp=1+n.*mcQ.*pbar.^(n-1).*(1-pbar).^2 - 2.*mcQ.*pbar.^2.*(1-pbar);
bdQdz=gamma./(1+D*M./(D+Qbar).^2);
bdcdz=-D*M./(D+Qbar).^2.*bdQdz;
ode = @(Qhat,chat) [1i.*w.*(X(1)./dQdp - X(2)) + X(2)./(D+pbar+cbar).*(1i.*w.*(D+pbar+cbar)-bdQdz)-X(1).*bdcdz, ...
X(2)./(D+Qbar+cbar).*(1i*w.*(D+pbar+cbar)-bdQdz)-X(1).*bdcdz];
ic = [bdQdz bdcdz+gamma];
[t,X] = ode45(ode,[0 1],ic);
function y = fphi(p,n,mcQ,Q)
y = p+mcQ*p^n*(1-p)^2-Q ;
end
Hey there!
I'm trying to solve the ODE above using the above IC's. Can anyone offer some help?
Adam Danz on 26 Feb 2021
@Brad Scott, in response to your flag, it would be better to improve the question by editing it rather than reposting it.

James Tursa on 23 Feb 2021
Edited: James Tursa on 23 Feb 2021
Just make your function handle return a column vector by using ; instead of , to separate the elements. E.g.,
ode = @(Qhat,X) [1i.*w.*(X(1)./dQdp - X(2)) + X(2)./(D+pbar+cbar).*(1i.*w.*(D+pbar+cbar)-bdQdz)-X(1).*bdcdz; ...
X(2)./(D+Qbar+cbar).*(1i*w.*(D+pbar+cbar)-bdQdz)-X(1).*bdcdz];
What is pbar?
Walter Roberson on 26 Feb 2021
My interpretation:
You are trying to solve a system of 1000 differential equations, but you only initialize with two boundary conditions.
If you passed in 2000 boundary conditions, interleaved, then at the end of the nested_odefun that I put in, put in
ode = ode(:);