# Subtracting results of NDSolve

I have 2 variables as a results from the NDSolve[] which are ua and uw. This is how I do it, and I got the result for ua and uw just fine and I can plot them.

eq1 = Derivative[2, 0][uw][z, t]*cwv -
Derivative[0, 1][ua][z, t]*cw == Derivative[0, 1][uw][z, t];

eq2 = Derivative[2, 0][ua][z, t]*cav -
Derivative[0, 1][uw][z, t]*ca == Derivative[0, 1][ua][z, t];

ic1 = {uw[zmax, t] == ui + p0, uw[z, tmin] == ui};

ic2 = {ua[z, tmin] == uatm};

sol2 = NDSolve[Flatten[{eq1, eq2, ic1, ic2}], {uw, ua}, {z, zmin, zmax}, {t, tmin,
tmax}, MaxSteps -> Infinity, MaxStepSize -> dt/step, InterpolationOrder -> All];


However, I want to subtract the two variables such that: result = (ua-101)-uw and I failed. I use this as a code to subtract

result = ua[z, t] /. sol2 - 101 - uw[z, t] /. sol2;

Can someone please tell me the reason why it is failed and what should I do to fix this?

The complete code is given as follows:

kw = 1*10^(-10);(*m/s*)

ka = 8*10^(-13);(*m/s*)

gw = 10;(*kN/m^3*)

m2w = 2.41*10^(-4);(*coefficient of compressibility of water phase with respect to matric suction*)

m1w = 0.8*10^(-4);(*coefficient of compressibility of water phase with respect to axial stress*)

m2a = 1.11*10^(-4);(*coefficient of compressibility with air phase with respect to matric suction*)

m1a = 0.37*10^(-4);(*coefficient of compressibility of air phase with respect to axial stress*)

p0 = 200;(*applied pressure -kPa*)

tmin = 0;(*minimum or time when the pressure is applied*)

tmax = 5000;(*latest inspected time*)

zmin = 0;(*top elevation*)

zmax = 0.02;(*bottom elevation*)

uatm = 101;(*kPa-atmospheric air pressure*)

step = 1000;(*The higher the step, the higher the accuracy*)

ui = -300;(*Initial matric suction*)

temp = 300;(*temperature in kelvin*)

wa = 28.97;(*kg/kmol*)

s = 0.7887;

n = 1.0696/(1 + 1.0696);

(*CALCULATION*)

r = 8.31432;(*Universal gas constant in j/(mol K*)

g = 9.8;(*gravity acceleration-m2/s*)

da = ka/g;

cwv = kw/(gw*m2w);

cw = ((1 - m2w/m1w)/(m2w/m1w));

ca = (m2a/m1a)/(1 - (m2a/m1a) - (1 - s)*n/(uatm*m1a));

cav = (da/(wa/(r*temp)*uatm))/(m1a*(1 - m2a/m1a) - (1 - s)*n);

dt = tmax - tmin;

(*Assuming that dua/dt is insignificant, Consolidation equation for unsaturated soil is duw/dt=cwv(d^2uw/dy^2)*)

eq = Derivative[2, 0][uw][z, t]*cwv == Derivative[0, 1][uw][z, t];

ic = {uw[zmax, t] == ui + p0, uw[z, tmin] == ui};

sol = NDSolve[Flatten[{eq, ic}], uw, {z, zmin, zmax}, {t, tmin, tmax},MaxSteps -> Infinity, MaxStepSize -> dt/step,InterpolationOrder -> All]

is = 300;(*ImageSize*)

Plot3D[uw[z, t] /. sol, {z, zmin, zmax}, {t, tmin, tmax}, PlotRange -> All, AxesLabel -> {z, t, uw[z, t]}, BoundaryStyle -> Thick, ImageSize -> 500]

Column[Manipulate[Plot[uw[z, t] /. sol, {z, zmin, zmax}, PlotRange -> All, AxesLabel -> {z, uw[z, t]}, ImageSize -> is], {t, tmin, tmax}], Manipulate[Plot[uw[z, t] /. sol, {t, tmin, tmax}, PlotRange -> All, AxesLabel -> {t, uw[z, t]}, ImageSize -> is], {z, zmin, zmax}],  Manipulate[Plot[{uw[z, tmin] /. sol, uw[z, t1] /. sol, uw[z, t2] /. sol, uw[z, t3] /. sol, uw[z, t4] /. sol, uw[z, tmax] /. sol}, {z, zmin,zmax}, AxesLabel -> {z, uw[z, t]}, ImageSize -> is], {t1, tmin + dt/300, tmin + dt/51}, {t2, tmin + dt/50, tmin + dt/16}, {t3, tmin + dt/15, tmin + dt/3}, {t4, tmin + dt/3, tmin + dt/1}]]

(*Assuming that dua/dt is significant*)
eq1 = Derivative[2, 0][uw][z, t]*cwv -
Derivative[0, 1][ua][z, t]*cw == Derivative[0, 1][uw][z, t];

eq2 = Derivative[2, 0][ua][z, t]*cav -
Derivative[0, 1][uw][z, t]*ca == Derivative[0, 1][ua][z, t];

ic1 = {uw[zmax, t] == ui + p0, uw[z, tmin] == ui};

ic2 = {ua[z, tmin] == uatm};

sol2 = NDSolve[Flatten[{eq1, eq2, ic1, ic2}], {uw, ua}, {z, zmin, zmax}, {t,tmin,tmax}, MaxSteps -> Infinity, MaxStepSize -> dt/step, InterpolationOrder -> All];

(*Coupled equation version*)
Column[Plot3D[uw[z, t] /. sol2, {z, zmin, zmax}, {t, tmin, tmax}, PlotRange -> All,AxesLabel -> {z, t, uw[z, t]}, BoundaryStyle -> Thick, ImageSize -> 500],
Plot3D[ua[z, t] /. sol2, {z, zmin, zmax}, {t, tmin, tmax}, PlotRange -> All, AxesLabel -> {z, t, ua[z, t]}, BoundaryStyle -> Thick, ImageSize -> 500]]

Column[Manipulate[Plot[uw[z, t] /. sol2, {z, zmin, zmax}, PlotRange -> Full, AxesLabel -> {z, uw[z, t]}], {t, tmin, tmax}],Manipulate[Plot[uw[z, t] /. sol2, {t, tmin, tmax},PlotRange -> Full, AxesLabel -> {t, uw[z, t]}], {z, zmin,zmax}]]Column[Manipulate[Plot[ua[z, t] /. sol2, {z, zmin, zmax}, PlotRange ->Full,AxesLabel -> {z, ua[z, t]}], {t, tmin, tmax}],Manipulate[Plot[ua[z, t] /. sol2, {t, tmin, tmax}, PlotRange -> Full, AxesLabel -> {t, ua[z, t]}], {z, zmin, zmax}]]

Column[Manipulate[Plot[{uw[z, tmin] /. sol2, uw[z, t1] /. sol2, uw[z, t2] /. sol2,uw[z, t3] /. sol2, uw[z, t4] /. sol2, uw[z, tmax] /. sol2}, {z, zmin, zmax}, AxesLabel -> {z, uw[z, t]}, PlotRange -> Full, AxesLabel -> {z, uw[z, t]}], {t1, tmin + dt/1000, tmin + dt/51}, {t2, tmin + dt/50, tmin + dt/16}, {t3, tmin + dt/15,tmin + dt/3}, {t4, tmin + dt/3, tmin + dt/1}], Manipulate[Plot[{ua[z, tmin] /. sol2, ua[z, t1] /. sol2, ua[z, t2] /. sol2, ua[z, t3] /. sol2, ua[z, t4] /. sol2, ua[z, tmax] /. sol2}, {z, zmin, zmax}, AxesLabel -> {z, ua[z, t]}, PlotRange -> Full], {t1, tmin + dt/500, tmin+ dt/51}, {t2, tmin + dt/50, tmin + dt/16}, {t3, tmin + dt/15, tmin + dt/3}, {t4, tmin + dt/3, tmin + dt/1}]]

(*Comparison*)

Manipulate[Plot[{uw[z, t] /. sol, uw[z, t] /. sol2}, {t, tmin, tmax},PlotRange -> All, AxesLabel -> {t, uw[z, t]}], {z, zmin, zmax}];

(*Matric Suction*)

result = ua[z, t] /. sol2 - 101 - uw[z, t] /. sol2;
Plot3D[result, {z, zmin, zmax}, {t, tmin, tmax}, PlotRange -> All, AxesLabel -> {z, t, \[Psi][z, t]}, BoundaryStyle -> Thick, ImageSize -> 500]


Regards,

-
This question might be useful for future visitors if you strip it down to essentials (in fact the bit @AlbertRetey points out for you). Please consider editing. – Yves Klett Feb 8 '13 at 7:50

result = (ua[z, t] /. sol2) - (uw[z, t] /. sol2)

as ReplaceAll (== /.) will apply the replacement rules at every level, you could just as well do:
result = ua[z, t] - uw[z, t] /. sol2