Question: Issue with pdsolve

Hi,

While using pdsolve for a coupled system of pdes, the maple results doen't match with the desired ones.

Here is the system

restart:

Eq1:=diff(f(eta,tau),eta$3)+(f(eta,tau)+h(eta,tau))*diff(f(eta,tau),eta$2)-diff(f(eta,tau),eta$1)^2

+(1+epsilon*cos(Pi*tau))*theta(eta,tau)=Omega*diff(f(eta,tau),tau);

Eq2:=diff(h(eta,tau),eta$3)+(f(eta,tau)+h(eta,tau))*diff(h(eta,tau),eta$2)-diff(h(eta,tau),eta$1)^2

+c*(1+epsilon*cos(Pi*tau))*theta(eta,tau)=diff(h(eta,tau),tau);

Eq3:=1/Pr*diff(theta(eta,tau),eta$2)+(f(eta,tau)+h(eta,tau))*diff(theta(eta,tau),eta$1)=Omega*diff(theta(eta,tau),tau);

ICs := {f(eta,0)=0,theta(eta,0)=1,h(eta,0)=0};

bc:={f(0,tau)=0,D[1](f)(0,tau)=0,h(0,tau)=0,D[1](h)(0,tau)=0,theta(0,tau)=1,D[1](f)(N,tau)=0,D[1](h)(N,tau)=0,theta(N,tau)=0};

I want to plot diff(theta(eta, tau), eta) at eta=0, for this I have rewritten the system as 

sys := {(diff(Nux(eta,tau), eta))/Pr+(f(eta, tau)+h(eta, tau))*Nux(eta, tau) = Omega*(diff(theta(eta, tau), tau)), diff(f(eta, tau), eta, eta, eta)+(f(eta, tau)+h(eta, tau))*(diff(f(eta, tau), eta, eta))-(diff(f(eta, tau), eta))^2+(1+epsilon*cos(Pi*tau))*theta(eta, tau) = Omega*(diff(f(eta, tau), tau)), diff(h(eta, tau), eta, eta, eta)+(f(eta, tau)+h(eta, tau))*(diff(h(eta, tau), eta, eta))-(diff(h(eta, tau), eta))^2+c*(1+epsilon*cos(Pi*tau))*theta(eta, tau) = diff(h(eta, tau), tau),diff(theta(eta, tau), eta)=Nux(eta,tau)};

Pr:=0.72:Omega:=0.2:c:=0.5:N:=5:epsilon:=0.0:pds:= pdsolve(sys,ICs union bc,numeric,spacestep=0.1,time=tau,timestep=0.01,range=0..N,errorest=true):p1:=pds:-plot(-Nux(eta,tau), eta=0, tau= 0..2):epsilon:=0.1:pds1:= pdsolve(sys,ICs union bc,numeric,spacestep=0.1,time=tau,timestep=0.01,range=0..N,errorest=true):p2:=pds1:-plot(-Nux(eta,tau), eta=0, tau= 0..2):epsilon:=0.3:pds2:= pdsolve(sys,ICs union bc,numeric,spacestep=0.1,time=tau,timestep=0.01,range=0..N,errorest=true):p3:=pds2:-plot(-Nux(eta,tau), eta=0, tau= 0..2):epsilon:=0.5:pds3:= pdsolve(sys,ICs union bc,numeric,spacestep=0.1,time=tau,timestep=0.01,range=0..N,errorest=true):p4:=pds3:-plot(-Nux(eta,tau), eta=0, tau= 0..2):plots[display]({p1,p2,p3,p4});

But the output doesn' match.

 pdeplot.mw

Thanks