@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.

@patient x0 (my x00) is to be treated as a parameter, so must be a name. Thus you should leave out the assignment to x0.
Leaving that out everything works.

@patient I had a look at your worksheet called exact_solution.mw.

I suggest shortening what you are doing there to:

restart;
u1 := -(-1-(1/4)*x^2/t)/(-t)^(3/2);
u2 := exp((1/8)*x^2/t);
F:=unapply(piecewise(x^2/t>-4,u1*u2,0),x,t);
plot([seq(F(x,-0.01*i),i=1..10)],x=-1..1,color=black);
##OK the last curve is not red, but in this animation it is:
plots:-animate(plot,[F(x,t),x=-1..1,color=`if`(t<-0.01,black,red)],t=-.1..-0.01,frames=10,trace=9);

So are you claiming that U(x) = F(x,t) solves the system SYS obtained above with some appropriate S(x)?
And what is that S(x) supposed to be?
Why is x^2/t >-4 appearing as a break?

########################
An exact solution can be found from earlier remarks. You confused me by renaming v(z) to s(z).
sol:= {u(z)=(1+z/4)*exp(z/8), s(z)=exp(z/4)};
odetest(sol,sys);
simplify(%) assuming z::real; #OK
##So a solution to SYS is
SOL:=subs(s(z)=S(x),u(z)=U(x),z=x^2/t,sol);
odetest(SOL,SYS);
simplify(%) assuming real;

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.

@Kitonum This is somewhat shorter:
simplify(temp) assuming positive;
radnormal(%);

@Christopher2222 Sorry about not catching the obvious mistake right away. Assuming that up is the same as up on your figure in the worksheet, then the altitude is given by h=R*sin(theta(t)) when theta is measured from the positive x-axis as you indicate in your figure.

The potential energy is then m*g*h + C, where the constant C can be chosen arbitrarily.

Had the angle theta been measured from the positive y-axis then the altitude would be given by h=R*cos(theta(t)) and the potential energy would again be m*g*h + C. Keeping the initial conditions then the events trigger would be theta(t)=Pi/2 (and theta(t)=-Pi/2, but that won't ever happen for your initial condition).

@patient Clearly the initial conditions for S and U must be different from the initial conditions for s and u. The change of variable from z to x obviously has an effect.
I won't change your (to my mind) unfortunate name x0 for the initial z-value (z0 appears to me to be more natural).
I shall call the x-value corresponding to z=x0 for x00.
That value is found from x^2=t*z to be x00 = sqrt(t*x0) if we choose x >0, which I shall do.

The derivative of S at x = x00 is found by using the chain rule to be
S'(x) = s'(z)*2*x/t
thus
S'(x00)= s'(x0)*2*x00/t.

Now we notice that the parameter t appears in the initial conditions as S(sqrt(t*x0)) = ...
It turns out that dsolve doesn't accept that, whereas the form S(x00) = ... with x00 as a parameter is accepted. Whether that should be considered a bug or not is irrelevant right now: we can work around it by using x00 as parameter instead of t.

Thus we can do:
SYS2:=subs(t=x00^2/x0,SYS);
ICS:={S(x00)=exp(x0^2-16),U(x00)=9*exp(x0^2-10),D(S)(x00)=-4*exp(-12)*2*x0/x00};
res:=dsolve(SYS2 union ICS, numeric, method = rkf45_dae, parameters = [x00]);
res(parameters=[sqrt(x0*(-.1))]); # This means that t = -0.1
plots:-odeplot(res,[[x,100*S(x)],[x,U(x)]],0..10);

Notice that I have multiplied S by 100 so it doesn't get drowned by U.

The procedure p can be modified as follows:
p := proc (t)
   res(parameters = [sqrt(t*x0)]);
   plots:-odeplot(res, [x, U(x)], 0 .. 4)
end proc;

plots:-animate(p, [t], t = -.1 .. -0.1e-1,view=0..0.05);

@Christopher2222 The zero for potential energy can be set arbitrarily. It is of no consequence for the equations of motion. In fact just look at the case in hand. L = T-U. You only use two derivatives of L. Both would eliminate an additive constant to U. The main thing is that since the gravitational pull is downwards, U should increase with the altitude.

@Christopher2222 Obviously the potential energy should increase with the altitude. Without a change of sign it doesn't. Whether U is positive or negative is irrelevant.