[求助]变系数微分方程求解

[356days]

问题是这样的：


俺编的程序是这样，请问各位大侠可行否？

function F=force(H,ro1,ro2,mu,tau,h0,Gv,lambda)
global rho1 rho2 mu tau g h0 Gv lambda;
F=81.2*(rho1*mu*tau^3)^(1/2)*a^3-4/3*pi*a^3*g*(rho2-rho1)+2*pi*h0^2*Gv*a/((H+1)*a)*(1+5.32*(H+1)*a/lambda)^(-1)

% F: the net force on particle. crossflow lift=81.2*(ro1*mu*tau^3)^(1/2)*a^3, Gravity=4/3*pi*a^3*g*(ro2-ro1), van der Waals=2*pi*h0^2*Gv*a/((H+1)*a)*(1+5.32*(H+1)*a/lambda)^(-1)
% rho1: the density of water
% rho2: the particle density
% mu: the solution viscosity
% tau: the wall shear rate
% g: the gravitation constant
% h0: the minimun interfacial separation distance
% Gv: van der Waals interfacial free energy
% lambda: characteristic decay length for van der Waals
% H: the dimentionless distance from particle to the surface of membrane. H=z/a-1 %(z is the vertical distance from particle to surface of membrane)
end

function f=f1(H)
f=sqrt(sin((H+1)*a))/(((quad('sqrt(sin((H+1)*a))',0,((H+1)*a))))^(1/3))

% f1(H): the universal function to account for hydrodynamice interaction according to Smoluchowski-Levich approximation.
end

function funb1=B1(H)
funb1=diff(log(f1(H)),H)-F(H)
end

function funb2=B2(H)
funb2=F(H)*(diff(log(f1(H))),H)-(diff(F(H),H))/f1(H)
end

function Con=Con-Frofile(H,y)
% the DE y'' + B1(H)y’ + B2(H)y = 0 written as a system.
% y is the concentration profile of particle dependent on H.
Con(1)=y(2);
Con(2)=-B1(H)*y(2)-B2(H)*y(1);
f=f(;
return

% trying to solve y''+ B1(H)y'+ B2(H)y = 0 at y(0)=0 y(2B/a-2)=0

clear all
c=1; % the first guess for y'(0)
Hu=2B/a-2; % the top of flow cell
global B; % the half hight of flow cell
tspan=[0 Hu];
while c > -999 % loop through different c values
c=input('enter c, -1000 to stop');
if c==-1000 break; end % if c=-1000 jump out of loop
y0=[0 c]; % initial condition [y(0) y'(0)]
[H,y]=ode45('Con-Frofile',tspan,y0);
plot(H,y(:,1)) % plot current solution
hold on % keep the plot
lgthy=length(y); % find the length of the vector y
disp([y(lgthy,1)]) % show y(H=Hu) value
cc=num2str(c); % convert value of c to a string variable
text(2.1,y(lgthy,1),cc) % write value of c on rhs of graph
end
text(2.1,5,'c guess') % write c guess on the plot at (2.1,5)
hold off
print -deps Con-Frofile