matlab中ode15s函数的用法,matlab ode15s方程中两个参数怎么调用其他函数带入求解...

本人用ode15s算一个偏微分方程组，该方程组中有两个参数要通过另外一个隐函数方程组求解，不知道怎么调用，ode15s求该解偏微分方程组代码([t,y]=ode15s(@fangcheng,tspan,y0,options);)如下：

function dydt=fangcheng(t,y)

kf=200;a=0.363;u=0.361032/a;c0=7.03;L=8.31e-2;N=20;dz=L/N;p=436.8/(1-0.363);

q1=y(1:N);c1=y(N+1:2*N);q2=y(2*N+1:3*N);c2=y(3*N+1:4*N);

%-----------定义dq/dt-----------------------------

for j=1:N

dq1dt(j)=kf/p*(c1(j)-ce1(j));

dq2dt(j)=kf/p*(c2(j)-ce2(j));

end

%---------定义dC/dt-----------------------------------

dc1dt(1)=-u*(c1(1)-c0)/dz-(1-a)/a*p*dq1dt(1);

dc2dt(1)=-u*(c2(1)-c0)/dz-(1-a)/a*p*dq2dt(1);

for j=2:N

dc1dt(j)=-u*(c1(j)-c1(j-1))/dz-(1-a)/a*p*dq1dt(j);

dc2dt(j)=-u*(c2(j)-c2(j-1))/dz-(1-a)/a*p*dq2dt(j);

end

dydt=[dq1dt dc1dt dq2dt dc2dt]';

其中ce1(j)、ce2(j)、q1(j)和q2(j)的关系是q1(j)=0.45967ce1(j)^0.3839/(ce1(j)^0.0457+0.3*ce2(j)^0.0003206，q2(j)=0.17356*ce2(j)^(-1.55201)/(ce2(j)^(-2)+0.1*ce1(j)^(-1))

如果用q1、q2表示c1、c2带入，程序是可以解出的，但是这是个隐函数方程组，不知道该怎么办，求大神相助