Problem with solving the inhomogeneous systems of ODE

Hi, I 'm trying to solve a system of ODE neumerically, in these equaions I have a funcion of time ,Pc, which was calculated seperately . I dont know how enter this function that now is a matrix of numbers to this system of odes. I know it s better to enter it with time dependent form but its form is complex and I m not able to extract its coeffisient.In fact so far I solve systems with constant inhomogeneous part .

Sorry I have a lot of constants which makes code long

I am new to MATLAB and would really appreciate the help.

With best wishes

Raha

function [t,y] = solveODEPcNew
%% constant
t=linspace(1e-11,-3,1000);
e=1.6*1e-19;
mili=1e-3;

nA=3.2;
nB=3.18;
nP=3.19;
c=3*1e8;
sigma=0.75;
%% structure
R1=24;%k/w
R2=7;
R3=9.4;
micro=1e-6;
lA=300*micro;
wA=micro;
tA=0.15*micro;%tickness
vA=lA*wA*tA;%volume
aA=lA*wA;%area
%pc structure
lP=250*micro;
wP=1.5*micro;
tP=0.25*micro;%tickness
vP=lP*wP*tP;
aP=lP*wP;
%Bragg section
lB=300*micro;
leffB=135*micro;
wB=1.5*micro;
tB=0.25*micro;%tickness
vB=lB*wB*tB;
aB=lB*wB;
RHOv=4.825e3;
V=[vA,vP,vB,vA,vP,vB];
Cv=3.124e2;
Cc=2.46e-7;
Cs=5.32e-5;
Cj=Cv*RHOv.*V;
%average Heat Capacity
C=mean(Cj);
T0=300;
% Ia=70*mili;
% Iast=70*mili;
% IB=20*mili;
% IBst=20*mili;
% Ip=2.7*mili;
% Ipst=2.7*mili;
%% ohmic heat term in sections Pi=a1i*Ii+a2i*Ii^2
rhoN=[9.24e-5,9.24e-5,2.19e-3,2.229e-3];
%thicness of differenr layer of 6 section
dna=[0.15,0.15,2.5,0.5];
dna=dna.*micro;
Ea=0.83;
% a1a=sigma*Ea;
% a2a=sigma^2*sum(rhoN.*dna)/aA;
% Rna=sum(rhoN.*dna)/aA;
%
dnB=[0.25,0.15,2.5,0.5];
dnB=micro.*dnB;
EB=0.95;
a1B=sigma*EB;
a2B=sigma^2*sum(rhoN.*dnB)/aB;
RnB=sum(rhoN.*dnB)/aB;
%
dnp=[0.25,0.15,2.5,0.5];
dnp=micro.*dnp;
Ep=0.95;%Ea=Ea*
l0=lA+lB+lP;%chip &substrate
rho0=2.229e-3;
d0=100*micro;
w0=400*micro;
rc=(rho0*d0)/(w0*l0);
a1a=sigma*Ea;
a2a=sigma^2*sum(rhoN.*dna)/aA;
IB=20*mili;
JB=sigma*IB/(e*vB);
Ba=-5.97*1e-27;
Ntr=1.37*1e24;
Na0=Ntr
NB0=Ntr;
Np0=Ntr;
Naini =3.1642e+24% with 70 mili
NBini=3e24;
Npini=4e24;
nAini=nA+Ba*(Naini-Na0);
nBini=nB+Ba*(NBini-NB0);
nPini=nP+Ba*(Npini-Np0);
Leffcini=nAini*lA+nBini*leffB+nPini*lP;
mNum= 2*Leffcini/(1553*1e-9);
LambdaB=1553*1e-9/(2*nBini);
g = 1.454401546848779e-13;
Sssini =1.7418e+22;
gammaN =3.8180e+08;
dalphadN=1.8*1e-21;%m^2
AprimA=zetaA*dalphadN;%A'a
k1=mNum*LambdaB/(Leffcini*g)*(-gammaPini*Ba+c*AprimA);
k2=mNum*LambdaB*lA/(Leffcini*lP)*(-gammaPini*Ba+c*AprimA)/g+(leffB-mNum*LambdaB)/lP;
Naini =3.1642e+24% with 70 mili
Nass = 3.8658e+24
tauB=1.0921e-08;% linear assumption
Tw=302;
b1=0.28*1e8;
b2=(0.28*1e-16-0.48*1e-17);
b3=(1+(1.6*1e-2)*(Tw-T0))*(0.52*1e-41);
%Rp(Np)
p1=b1;
p2=b2;
p3=b3;
b3=(1+(1.6*1e-2)*(Tw-T0))*(0.52*1e-41);
JBini=b1*NBini+b2*NBini^2+b3*NBini^3;
Jpini=p1*Npini+p2*Npini^2+p3*Npini^3;
Iaini=70*mili;
%Jaini=sigma/e*vA*Iaini
LambdaB=1553*1e-9/(2*nBini);
%Na=Naini+k1*tauB*(JB-JBini)*(1-exp(-t./tauB));
%Ja=Ndota+gammaN.*Na+g.*(Na-Ntr).*Sssini;
Jp=Jpini+k2*(JB-JBini);
k1=mNum*LambdaB/(Leffcini*g)*(-gammaPini*Ba+c*AprimA);
%k2=mNum*LambdaB*lA/(Leffcini*lP)*(-gammaPini*Ba+c*AprimA)/g+(leffB-mNum*LambdaB)/lP;
Na=Naini+k1*tauB*(JB-JBini)*(1-exp(-t./tauB));
Ndota=k1*(JB-JBini)*exp(-t./tauB);
Ja=Ndota+gammaN.*Na+g.*(Na-Na0).*Sssini;
Ia=e*vP.*Ja/sigma;
%
Iamax=0.212;
Ipmax=Ip;
IBmax=IB;
PaBar=a1a*Iamax+a2a*Iamax^2;
Pa=a1a.*Ia+a2a.*Ia.^2;
Iast=(-a1a+(a1a^2+sqrt(a1a^2-4*a2a.*(Pa-2*PaBar))))/(2*a1a);
% Past=a1a*Iast+a2a*Iast^2;
IBst=20*mili;
Ip=2.7*mili;
Ipst=2.7*mili;
I=Ia+Iast+IB+IBst+Ip+Ipst;
%P=Pa+Past+PB+PBst+Pp+Ppst;
Ileak=(1-sigma).*I.^2;
Pc=rc*Ileak.^2;
%% Matrix of coeffeisients

m11=-1/(R2*Cc);
m12=-m11;
m13=1/(Cc*R1);
% m14=0;
m21=1/(R2*Cs);
m22=-(1/R2+1/R3)*(1/Cs);
m23=0;
% m24=0;
m31=6/(R2*Cc);
m32=-6/(R2*Cc);
m33=-6/(Cc*R1)-1/(R1*C);
%% inhomogeneous part
U1=Pc/Cc;
U2=0;
U3=-6*Pc/Cc;%
%% solve the system T0
tspan=[1e-11 1e-3];
y0=[0 0 0 ];
[t,y] = ode15s(@odefcn,tspan,y0);
%% ODE function
    function dydt =  odefcn(~,y)
        dydt=zeros(3,1);
        dydt(1)= m11*y(1)+ m12*y(2)+ m13*y(3)+U1;
        dydt(2)= m21*y(1)+ m22*y(2)+ m23*y(3)+U2;
        dydt(3)= m31*y(1)+ m32*y(2)+ m33*y(3)+U3;
    end
end

[Alan Stevens]

You have some undefined parameters: zetaA, gammaPini and Ip.

Your large matrices such as I, Ia etc. seem to be mainly NaNs or infs.

[Walter Roberson]

Your t is a vector of values, which forces a number of other variables to be a vector of values, including U1 and U3. So in

        dydt(1)= m11*y(1)+ m12*y(2)+ m13*y(3)+U1;
        dydt(2)= m21*y(1)+ m22*y(2)+ m23*y(3)+U2;
        dydt(3)= m31*y(1)+ m32*y(2)+ m33*y(3)+U3;

the first and the third of those wants to be a vector of values, one for each different t; meanwhile the second of the lines only wants to be a scalar.

As your U1 and U3 are built from "time", perhaps you should be using the original t vector together with the value passed in as the first parameter to odefcn to interpolate a particular U1 and U3 value ?

[raha ahmadi]

Dear Walter

Thanks for your help. I realized my problem and tried to change it . I did a long calculation and replace U1 with its equation (in fact I replace U1 and U3 with their equations as a function of time.after that I introduced time ,t, as a 4th input as you suggested.

I really appreciated if you check the new code. Also now I face the other error

"Warning: Failure at t=1.000000e-07. Unable to meet integration tolerances without

reducing the step size below the smallest value allowed (2.117582e-22) at time t.

> In ode15s (line 730)

In solveODEPcNew (line 133) " this means ode15s cant solve it.

Thanks in advance

Raha

Dear Alan Stevens

sorry I fotgot send the other parts of my code. my main problem is what Walter Roberson said. thanks for time

function [t,y] = solveODEPcNew
%% constant
mili=1e-3;
tauB=1.0921e-08;
AprimA =5.850000000000000e-22;
sigma=0.75;
R1=24;%k/w
R2=7;
R3=9.4;
micro=1e-6;
lA=300*micro;
wA=micro;
tA=0.15*micro;%tickness
vA=lA*wA*tA;%volume
aA=lA*wA;%area
lP=250*micro;
wP=1.5*micro;
tP=0.25*micro;%tickness
vP=lP*wP*tP;
aP=lP*wP;
lB=300*micro;
wB=1.5*micro;
tB=0.25*micro;%tickness
vB=lB*wB*tB;
aB=lB*wB;
V=[vA,vP,vB,vA,vP,vB];
Cv=3.124e2;
Cc=2.46e-7;
Cs=5.32e-5;
RHOv=4.825e3;
Cj=Cv*RHOv.*V;
C=mean(Cj);
Ntr=1.37*1e24;
Na0=Ntr
l0=lA+lB+lP;%chip &substrate
rho0=2.229e-3;
d0=100*micro;
w0=400*micro;
rc=(rho0*d0)/(w0*l0);
x=rc*(1-sigma)^2
Ba=-5.97*1e-27;
c=3e8;
IB=20*mili;
IBst=20*mili;
Ip=2.7*mili;
Ipst=2.7*mili;
Ires=IB+IBst+Ip+Ipst;
Ea=0.83;
a1a=sigma*Ea;
rhoN=[9.24e-5,9.24e-5,2.19e-3,2.229e-3];
dna=[0.15,0.15,2.5,0.5];
dna=dna.*micro;
a2a=sigma^2*sum(rhoN.*dna)/aA;
lambda0=1553e-9
dndomega=1e-16;
nA=3.2;
freq=c/(nA*lambda0);
omega=2*pi*freq;
zetaA=0.325;
zetaB=0.488;
zetaP=0.488;
ngA=nA+dndomega*omega;
dgdNa=4.83*1e-21;
VgA=c/ngA
g=zetaA*dgdNa*VgA
Sssini = 1.648956271484272e+22;
gammaN =3.8180e+08;
mNum =2.805341422150676e+03
Naini =3.1642e+24
tauB=1.0921e-08;
JB =8.333333333333332e+32;
%JBini=b1*NBini+b2*NBini^2+b3*NBini^3;
JBini =4.376928000000000e+32;
DeltaJB=JB-JBini;
Leffcini = 0.002178347614300;
gammaPini =3.339100000000000e+11;
NBini=3e24;
nB=3.18;
Npini=4e24;
NB0=Ntr;
nBini=nB+Ba*(NBini-NB0);
LambdaB=1553*1e-9/(2*nBini);
k1=mNum*LambdaB/(Leffcini*g)*(-gammaPini*Ba+c*AprimA);
T0=300;
Iamax=0.212;
Ipmax=Ip;
IBmax=IB;
PaBar=a1a*Iamax+a2a*Iamax^2;
EB=0.95;
a1B=sigma*EB;
dnB=[0.25,0.15,2.5,0.5];
dnB=micro.*dnB;

a2B=sigma^2*sum(rhoN.*dnB)/aB;
PB=a1B*IB+a2B*IB^2;
Ep=0.95;
a1p=sigma*Ep;
dnp=[0.25,0.15,2.5,0.5];
dnp=micro.*dnp;
a2p=sigma^2*sum(rhoN.*dnp)/aP;
Pp=a1p*Ip+a2p*Ip^2;
Ppst=a1p*Ipst+a2p*Ipst^2;
PBst=PB;
P=2*PaBar+PB+PBst+Pp+Ppst;
d1=DeltaJB*k1*(1-tauB*(gammaN+g*Sssini));
d2=(gammaN+g*Sssini)*(Naini+k1*tauB*DeltaJB)-g*Sssini*Na0;
d3=a1a*d1+2*a2a*d1*d2;
d4=a2a*d1^2;
d5=a1a*d2+a2a*d2;
y01=(R2+R3)*P+T0;
y02=R3*P+T0;
y03=R1*P;
%%
m11=-1/(R2*Cc);
m12=-m11;
m13=1/(Cc*R1);
% m14=0;
m21=1/(R2*Cs);
m22=-(1/R2+1/R3)*(1/Cs);
m23=0;
% m24=0;
m31=6/(R2*Cc);
m32=-6/(R2*Cc);
m33=-6/(Cc*R1)-1/(R1*C);

% U1=1/Cc*(x*(d1*exp(-t/tauB)+d2-a1a/(2*a2a)+1/(2*a2a)*sqrt(a1a^2-4*a2a*d5+...
%     8*a2a*PaBar-4*a2a*d3*exp(-t/tauB)-4*a2a*d4*exp(-2*t/tauB))+Ires)^2)%Pc/Cc;
% U3=-6*1/Cc*(x*(d1*exp(-t/tauB)+d2-a1a/(2*a2a)+1/(2*a2a)*sqrt(a1a^2-4*a2a*d5+...
%     8*a2a*PaBar-4*a2a*d3*exp(-t/tauB)-4*a2a*d4*exp(-2*t/tauB))+Ires)^2);%
%% solve the system T0
tspan=[1e-11 1e-3];
y0=[y01 y02 y03 1e-11];
[t,y] = ode15s(@odefcn,tspan,y0);
%% ODE function
    function dydt =  odefcn(~,y)
        dydt=zeros(4,1);
        dydt(1)= m11*y(1)+ m12*y(2)+ m13*y(3)+1/Cc*(x*(d1*exp(-y(4)/tauB)+d2-a1a/(2*a2a)+1/(2*a2a)*sqrt(a1a^2-4*a2a*d5+...
    8*a2a*PaBar-4*a2a*d3*exp(-y(4)/tauB)-4*a2a*d4*exp(-2*y(4)/tauB))+Ires)^2);
        dydt(2)= m21*y(1)+ m22*y(2)+ m23*y(3);
        dydt(3)= m31*y(1)+ m32*y(2)+ m33*y(3)+-6*1/Cc*(x*(d1*exp(-y(4)/tauB)+d2-a1a/(2*a2a)+1/(2*a2a)*sqrt(a1a^2-4*a2a*d5+...
    8*a2a*PaBar-4*a2a*d3*exp(-y(4)/tauB)-4*a2a*d4*exp(-2*y(4)/tauB))+Ires)^2);
    end
end