@javaneh0 You say that you forgot the initial condition "ph+pc+pa+pgd+prd+prj+pf=1".

Since you are setting
ics := ph(0) = 1, pc(0) = 0, pa(0) = 0, prj(0) = 0, prd(0) = 0, pgd(0) = 0, pf(0) = 0;

this condition is not additional, it is automatically satisfied.

Maybe you mean the requirement ph(t)+pc(t)+pa(t)+pgd(t)+prd(t)+prj(t)+pf(t)=1  for all t?

@javaneh0 My code above is plain text and can be copied and pasted into a fresh Maple worksheet without any problem.

If you have problems with that let me know.

@patient Sorry about venting my frustration with the document mode and 2D-math. It is not your fault. It is the default environment in Maple these days if you don't change it, which, however, is easily done under Tools/Options.

For the exact solution just use piecewise:
solp:={u(z)=piecewise(z>-4,(1+z/4)*exp(z/8),0),v(z)=piecewise(z>-4,exp(z/4), -4*exp(-1)*1/z)};
SOLp := subs(u(z) = (-t)^(3/2)*U(x), v(z) = -V(x)*t, z = x^2/t, solp);

SOL2p:=solve(SOLp,{U(x),V(x)});
odetest(SOL2p,SYS);
simplify(%) assuming x^2/t>-4;
simplify(%%) assuming x^2/t<-4;

plots:-animate(plot,[subs(SOL2p,[U(x),V(x)]),x=0..2,color=[red,green],thickness=3],t=-1..-0.1);
For the numerical solution I see no reason not to just use the exact solution after integration stops.
After all it is a somewhat artificial stop that is made. The exact solution for U is continuous, yes, but not differentiable at the point where it becomes zero.

resE:=dsolve(SYS2 union ICS,numeric,parameters=[x0],method=rkf45_dae,events=[[x=2*x0-.001,halt]]);
pE := proc (t, itvl::range) local p1,p2;
  resE(parameters = [sqrt(-t)]);
  p1:=plots:-odeplot(resE, [[x, U(x)], [x, V(x)]], itvl, _rest);
  p2:=plot([0,4*exp(-1)/x^2],x=2*sqrt(-t)..rhs(itvl),_rest);
  plots:-display(p1,p2)
end proc;
plots:-animate(pE,[t,0.001..2,thickness=3],t=-1..-.1);

@emmantop Expanding exp(beta/theta) in a taylor series raises the question:
With respect to what variable and about what point?

Also as I was saying above even replacing exp(beta/theta) by 1 doesn't help.

Even trying replacing exp(...) by 1 as in

dsolve(eval(sys, {epsilon = E[1],exp=1}),numeric, method = bvp[midrich],maxmesh=256);

doesn't work.

When solving for the fourth derivative of f and expanding there is a term having
f''(eta)*theta'(eta) / ( f(eta)*theta(eta)^2)
and the exp-factor and a factor beta.

I tried setting beta=0 instead of replacing exp directly by 1. This has the effect of removing that term besides having the desirable side effect of setting exp=1.

Then there was no problem:

P1 := {Ec = .5, N = 10, Re = 5,lambda = -1, n = 2, sigma = 10, k[1] = .9}; #No beta value
#ode1, ode2, and bcs as before.
sys1:=eval([ode1, ode2, bcs], P1);
res:=dsolve(eval(sys1, {epsilon = E[1],beta=0}),numeric, method = bvp[midrich],maxmesh=256);
plots:-odeplot(res,[[eta,f(eta)],[eta,theta(eta)]],0..10);

####Looking at the term left out by setting beta=0 (for what it is worth):
T:=diff(theta(eta),eta)/theta(eta)^2;
F:=diff(f(eta),eta,eta)/f(eta);
plots:-odeplot(res,[[eta,F*T]],0..10)
#Huge at the end of the interval. 
 

@patient I detest document mode and 2D math. I have to fight it wasting valuable time. Personally I never use it.

Here is the worksheet version in 1D-math.

I added a dsolve-version at the end using events to stop integration a small amount before x = 2*x0 at which the exact solution for U(x) will become zero. Since U(x) appears in the denominator for V''(x) =  ... this will be critical for a numerical solution. Furthermore plotting is done there from x=0.001 in order to avoid the problem at x=0.

MaplePrimes15-05-09odesys_parametersC.mw

I have no solution, but notice that if you remove the exp-terms then there is no problem.

As in

dsolve(eval(sys, {epsilon = E[1],exp=0}),numeric, method = bvp[midrich],maxmesh=256);

@patient I'm afraid we have returned to one of the problems I pointed out in my original answer: You cannot give initial conditions at z=0, which here corresponds to x=0.

Another unrelated point:
To check the exact solution given somewhat implicitly by SOL you need tro isolate U(x) and V(x):

solve(SOL,{U(x),V(x)});
odetest(%,SYS);
SOL2:=simplify(%) assuming real;

##############################
I now choose initial conditions taken from the exact solution.
I use z=-1 (arbitrarily: We need z<>0).
eval(SOL2,{x=x0,t=-x0^2});
ICS1:=simplify(%) assuming x0>0;
convert(diff(SOL2[2],x),D);
eval(%,x=x0);
ICS2:=eval(%,{t=-x0^2});
ICS:=ICS1 union {ICS2};
SYS2:=simplify(subs(t=-x0^2,SYS)) assuming x0>0;
numer(lhs(SYS2[2]))/x0^5;

res:=dsolve(SYS2 union ICS,numeric,parameters=[x0],method=rkf45_dae,maxfun=100000);
p := proc (t, itvl::range) res(parameters = [sqrt(-t)]); plots:-odeplot(res, [[x, U(x)], [x, V(x)]], itvl, _rest) end proc;
plots:-animate(p,[t,0..2,thickness=3],t=-1..-.1); p1:=%:
plots:-animate(plot,[subs(SOL2,[U(x),V(x)]),x=0..2,color=[red,green]],t=-1..-0.1); p2:=%:
plots:-display(p1,p2);

@javaneh0 After making the corrections that are pointed out by Makiyan then using method=laplace will give you a shorter looking result:

sys_ode := diff(ph(t),t) = (1-yc)*pc(t)+yh*prj(t)+urd*prd(t)+ugd*pgd(t)-yc*ph(t), diff(pc(t),t) = yc*ph(t)-(2-yc)*pc(t), diff(pa(t),t) = ya*pc(t)-pf(t), diff(prj(t),t) = yrj*pa(t)+prj(t), diff(prd(t),t) = -urd*prd(t)+yrd*pa(t), diff(pgd(t),t) = -ugd*pgd(t)+ygd*pa(t), diff(pf(t),t) = (1-ygd-yrj-yrd)*pa(t)+(1-yh)*prj(t);

ics := ph(0) = 1, pc(0) = 0, pa(0) = 0, prj(0) = 0, prd(0) = 0, pgd(0) = 0, pf(0) = 0;

vars:=indets({sys_ode},function(identical(t)));
dsolve({sys_ode,ics},vars,method=laplace);

To get a shorter or more concrete result you need to give concrete values to the parameters, i.e. to the names yc, yh, etc.
You may in that case consider solving numerically.
Example of solving numerically:
All parameters are set to the same value 1 (easy to do that is why):

params:=indets({sys_ode},name) minus {t}=~1;
dsolve({sys_ode,ics},vars,method=laplace);
res:=dsolve(eval({sys_ode,ics},params),numeric);
plots:-odeplot(res,[t,ph(t)],0..5);

@mobiusinfi You can only have 2 boundary conditions. So think again.
I tried:

restart;
k :=21.37595739:
c :=1.760652758:

v:=unapply(u(x,t)+k*u(x,t)^(2),x,t);

PDE := diff(u(x, t), t)-(diff(v(x, t), x, x)) = 0;
#Leaving one of your boundary conditions out:
IBC1 := {u(x, 0) = 0, D[1](v)(0, t) = 0, D[1](v)(1, t) = c};
S1 := pdsolve(PDE, IBC1, numeric, time = t);
S1:-animate(t=5);
#Leaving another out instead produces the error known from earlier:
IBC2 := {u(1,t)=1,u(x, 0) = 0, D[1](v)(0, t) = 0};
S2 := pdsolve(PDE, IBC2, numeric, time = t);
S2:-animate(t=5);

@tomleslie I doubt that the LambertW solution is going to solve the problem. I tried the following:


restart;
PDE := diff(u(x, t), t)-diff(u(x, t)+k*u(x, t)^2, x, x) = 0;
sol:=pdsolve(PDE);
#This solution appears to be of the form
sol2:=u(x,t)=A*LambertW(-exp(a*x+b*t+c))+B;
# where A, a,b,c, B are 5 arbitrary constants (not necessarily real!)
#However, their number is reduced to 3:
pdetest(sol2,PDE);
collect(numer(%),LambertW);
coeffs(subs(LambertW(-exp(a*x+b*t+c))=L,%),L);
solve({%},{B,b});
expand(%);
sol3:=eval(sol2,%);
#Clearly sol3 is not going to solve the initial value problem:
eval(sol3,t=0);

What we are getting is more likely just a building block as when pdsolve gives us in case k=0:
pdsolve(eval(PDE,k=0),build);

I looked at your pde in the following form, where your k = 21.37595739 and your c = 1.760652758.
Even using k=0 gave problems!

restart;
PDE := diff(u(x, t), t)-diff(u(x, t)+k*u(x, t)^2, x, x) = 0;
IBC := {u(1, t) = 1, u(x, 0) = 0, (D[1](u))(1, t) = c};

res:=pdsolve(eval(PDE,k=0.00001),eval(IBC,c=1),numeric,time=t,range=0..1);
res:-animate(t=.3);

What kind of behavior are you expecting? Your problem obviously comes from some physical situation. What is that?

Example:

You need the independent variable x in the ode

diff(y(x),x) = y(x)

You cannot just write

diff(y(x),x) = y

That is very likely your problem. There may be more problems.

@patient But your sol is not a solution of sys in your latest worksheet numeric_solution2.mw:

odetest(sol,sys);
simplify(%) assuming z::real,t<0;

You see that the second member is not identically zero.


Now try with coefficients A and B:
sol1 := {u(z) = (1+(1/4)*z)*exp((1/8)*z)*A, v(z) = exp((1/4)*z)*B};
odetest(sol1,sys);
simplify(%) assuming real;
collect(%,z);

Notice that we need A=sqrt(B)  (unless B=0).

So the question is: What do you want B to be?

@patient 
I wrote at the end of my reply Exact solution:
"The initial conditions which this solution satisfies at x00=sqrt(x0*t) are easily computed and there is a discrepancy between those and the ones I named ICS in your numerical approach."

Taking a look at that in the present context:

ICS; #Initial conditions you use for U and V
SOL; #The exact solution for U and V
#We find the initial conditions actually satisfied at x=x0=sqrt(t*z0) by the exact solution:
ICS0 := eval(SOL, x = sqrt(t*z0));
diff(SOL[2], x);
convert(%, D);
eval(%, x = sqrt(t*z0));
ICS1 := simplify(%) assuming t<0;
ICSE :=ICS0 union {ICS1};


You see that there is quite a discrepancy between ICS and ICSE.