# NAME    :- IBRAR AHMAD 
TUTORIAL:- GROUP 5TH   
ROLL NO :- 11232225    
DATE    :- 2-4-2013   

# 1. Lake Pollution Model

kill(all);
depends([S,Cin],t)$
eqn:'diff(S,t)=10*Cin-S/10;
gsoln:ode2(eqn,S,t);
ic1(gsoln,t=0,S=S0);

kill(all);
depends([S,Cin],t)$
eqn:'diff(S,t)=10*Cin-S/10;
gsoln:ode2(eqn,S,t);
psoln:ic1(gsoln,t=0,S=S0);
ev(subst([Cin=1,S0=1],psoln));
define(f(t),rhs(%o5));
wxplot2d(f(t),[t,0,15],[ylabel,"f(t)"],
[style,[lines,2]]);

kill(all);
depends([S,Cin],t)$
eqn:'diff(S,t)=10*Cin-S/10;
gsoln:ode2(eqn,S,t);
psoln:ic1(gsoln,t=0,S=S0);
ev(subst([Cin=0.3+0.1*sin(t),S0=1],psoln));
define(f(t),rhs(%o5));
wxplot2d(f(t),[t,0,15],[ylabel,"f(t)"],
[style,[lines,2]]);

kill(all)$
de:'diff(C,t)=(F/V)*c-(F/V)*C;
gsoln:ode2(de,C,t);
psoln:ic1(gsoln,C=C(0),t=0); /*C(0) is the initial concentration */
subst(c=0,psoln);
solve(rhs(%)=(5/100)*C(0),t);

kill(all)$
c:3;
V:28;
F:48;
de:'diff(C,t)=(F/V)*c-(F/V)*C;
gsoln:ode2(de,C,t);
psoln:ic1(gsoln,C=b,t=0);
for i:0 thru 6 step 1 do (g[i]:ev(rhs(psoln),b=i));
wxplot2d([g[0],g[1],g[2],g[3],g[4],g[5],g[6]],[t,0,8],[y,0,6],
[legend,"C=0","C=1","C=2","C=3","C=4","C=5","C=6"],
[xlabel,"t"],[ylabel,"c(t)10^6"]);

kill(all)$
ratprint:false$
c:4;
F:12*(1+6*sin(2*%pi*t));
V:28;
de:'diff(C,t)=(F/V)*c-(F/V)*C;
gsoln:ode2(de,C,t);
psoln:ic1(gsoln,C=b,t=0);
for i:0 thru 6 step 1 do (g[i]:ev(rhs(psoln),b=i));
wxplot2d([g[0],g[1],g[2],g[3],g[4],g[5],g[6]],[t,0,12],[y,0,8],
[legend,"C=0","C=1","C=2","C=3","C=4","C=5","C=6"],
[xlabel,"t"],[ylabel,"c(t)10^6"]);

kill(all)$
ratprint:false$
c:10+10*cos(2*%pi*t);
F:12;
V:28;
de:'diff(C,t)=(F/V)*c-(F/V)*C;
gsoln:ode2(de,C,t);
psoln:ic1(gsoln,C=b,t=0);
for i:0 thru 6 step 1 do (g[i]:ev(rhs(psoln),b=i));
wxplot2d([g[0],g[1],g[2],g[3],g[4],g[5],g[6]],[t,0,12],[y,0,12],
[legend,"C=0","C=1","C=2","C=3","C=4","C=5","C=6"],
[xlabel,"t"],[ylabel,"c(t)10^6"]);

kill(all)$
ratprint:false$
c:10+10*cos(2*%pi*t);
F:12*(1+6*sin(2*%pi*t));
V:28;
de:'diff(C,t)=(F/V)*c-(F/V)*C;
points:rk(rhs(de),C,0,[t,0,8,0.1])$
S0:makelist([points[i][1],points[i][2]],i,1,length(points))$
points:rk(rhs(de),C,1,[t,0,8,0.1])$
S1:makelist([points[i][1],points[i][2]],i,1,length(points))$
points:rk(rhs(de),C,2,[t,0,8,0.1])$
S2:makelist([points[i][1],points[i][2]],i,1,length(points))$
points:rk(rhs(de),C,3,[t,0,8,0.1])$
S3:makelist([points[i][1],points[i][2]],i,1,length(points))$
points:rk(rhs(de),C,4,[t,0,8,0.1])$
S4:makelist([points[i][1],points[i][2]],i,1,length(points))$
points:rk(rhs(de),C,5,[t,0,8,0.1])$
S5:makelist([points[i][1],points[i][2]],i,1,length(points))$
points:rk(rhs(de),C,6,[t,0,8,0.1])$
S6:makelist([points[i][1],points[i][2]],i,1,length(points))$
wxplot2d([[discrete,S0],[discrete,S1],[discrete,S2],[discrete,S3],
[discrete,S4],[discrete,S5],[discrete,S6]]);

# 2. Drug Assimliation in the blood

# Case I : Case of a Single cold pills:

kill(all)$
depends([x,y],t);
eqn1:'diff(x,t)=-k1*x;
eqn2:'diff(y,t)=k1*x-k2*y;
solx:ic1(ode2(eqn1,x,t),t=0,x=x0)$
g1:ev(solx,k1=1.3680,x0=1)$
g2:ev(solx,k1=0.6931,x0=1)$
wxplot2d([rhs(g1),rhs(g2)],[t,0,10],[legend,"Decongestant","Antihistamine"],
[style,[lines,1,1],[lines,1,2]],[ylabel,"x/x0"]);
f:ev(eqn2,solx)$
soly:ic1(ode2(f,y,t),t=0,y=0)$
g3:ev(soly,k1=1.3680,k2=0.1386,x0=1)$
g4:ev(soly,k1=0.6931,k2=0.0231,x0=1)$
wxplot2d([rhs(g3),rhs(g4)],[t,0,20],[legend,"Decongestant","Antihistamine"],
[style,[lines,1,1],[lines,1,2]],[ylabel,"y/y0"]);

# Case II : A course of Cold Pills:

kil(all)$
ratprint:false$
depends([x,y],t);
I:1$
eqn1:'diff(x,t)=I-k1*x;
eqn2:'diff(y,t)=k1*x-k2*y;
solx:ic1(ode2(eqn1,x,t),t=0,x=0)$
g1:ev(solx,k1=1.3860)$
g2:ev(solx,k1=0.6931)$
print("curve of medicine in GI tract is")$
wxplot2d([rhs(g1),rhs(g2)],[t,0,60],
[legend,"decongestant","Antihistamine"],
[style,[lines,2,1],[lines,2,2]],[xlabel,"t(hours)"],
[ylabel,"x/x0"]);
f:ev(eqn2,solx)$
f1:ev(f,k1=1.3860,k2=0.1386)$
f2:ev(f,k1=0.6931,k2=0.0231)$
f3:ic1(ode2(f1,y,t),t=0,y=0)$
f4:ic1(ode2(f2,y,t),t=0,y=0)$
print("curve of medicine  in the  bloodstream is")$
wxplot2d([rhs(f3),rhs(f4)],[t,0,60],
[legend,"decongestant","Antihistamine"],
[style,[lines,2,1],[lines,2,2]],[xlabel,"t(hours)"],
[ylabel,"y/y0"]);