The phase space of a 5th order ODE is five-dimensional.  The current version of Maple does three-dimensional pictures quite well.  You will have to wait for a future release of Maple to do plots in five dimensions.

@vercammen If you know about Maple's procs, then a better way would be through a proc, as in

make_prob := proc({n:=10, alpha:=-60,  beta:=20, P_val:=4})
    local Qs := alpha + beta*P;
    printf("There are %a identical firms.  "
           "Supply is given by Qs = %a.  "
           "When P = %a then quantity supplied is equal to %a.  "
           "Calculate the ...\n", n, Qs, P_val, eval(Qs, P=P_val));
end proc:

Then you may call the proc any number of times with the desired arguments.  Missing arguments take on default values as indicated in the first line of the proc's definition.  For example:

make_prob();   # all default args
There are 10 identical firms.  Supply is given by Qs = -60+20*P.  When P = 4 then quantity supplied is equal to 20.  Calculate the ...

make_prob(n=40, P_val=5);   # some optional args
There are 40 identical firms.  Supply is given by Qs = -60+20*P.  When P = 5 then quantity supplied is equal to 40.  Calculate the ...

@Neel There is no general way for that.  Each case needs to be treated on its own, but the steps are pretty much always the same..  So if you see what I did for the clamped-clamped case, you may do the same to analyze the clamped-free case (a cantilever beam).

Here is a pointer to simplify your calculations just a bit.

In what I presented earlier, I went along with general solution of the ODE produced by Maple as a linear combination of exponential and trigonometric functions with coefficients _C1 through _C4.  For our purposes that's not the best form of the general solution.  The first two terms may be replaced with hyperbolic functions, as in _C1*sinh(mu*x/L) + C_2*cosh(mu*x/L).  Since _C1 and _C2 are arbitrary, the changed expression is equivalent to the original.  I suggest that you make that change in the worksheet (manually) and then continue the calculations.  You will see that things work out more smoothly that way.

@Aminaple The idea suggested by Dr. Venkat Subramanian is quite good and I like it better than the solution that I presented earlier.  Here are the details of this alternative approach.

Download mw2.mw

@Neel 

@Neel What are the boundary condtions?  I cannot tell by reading through your worksheet.

@janhardo Your interpretation of my calculations is quite correct.  Well done!

@AHSAN Even after changing q(y)=j to u(y)=j the problem still makes no sense.  The "bc" and "bc1" in your formulation stand for "boundary condition" but you have not said where the boundaries are.

Go to where you found the problem, read it carefully, then come back with the correct statement.

I looked through your worksheets but was unable to understand them.  It will be very helpful to your readers to have some verbal description of the purpose of the worksheets.  You may want to supply the following information:

1. It appears that you are solving the PDE 
   
   with the initial conditions u=(x,0)=sin(Pi*x), diff(u(x,t),t)=0 on the interval (0.1).  Is that correct?
2. You state the boundary conditions as
   bc1:=u(x,t):        bc2:=u(x,t):
   I don't know that those mean.  Explain.
3. It appears that you discretize the equation through centered finite differencing in the x direction. Is that correct?
4. Then you discretize the equations in the time direction in some way.  That does not seem to be consistent with your stated objective of applying the method of lines, which approaches the problem through solving a system of ODEs, but perhaps I have misunderstood.  Explain.
5. What is the idea behind your non-classical method?
6. You write "I need to compare the result of the method with the exact solution".  What is the exact solution?
7. You must know that you may solve your PDE with Maple's pdsolve() (numerically).  You are aware of that, aren't you?
8. Why have you posted two worksheets?  What is the purpose of each?  Provide a guide.

@janhardo There are neither double integrals nor polar coordinates there.  The area is calculated by simple integration in Cartesian coordinates.

@AHSAN Your bc says q(y) = 0.  Your bc1 says q(y) = j .  Do you want both?

@Fancypants Yes, equation (20) looks good.

@Fancypants Okay, that's better.  Now you have an ODE for S[y](t).  The ODE involves phi(t), psi(t), theta(t), f(t), etc., with various subscripts.  Are those known/given functions?  If so, then your equation (13) is a simple second order linear ODE which is trivial to solve.  On the other hand, if you expect to determine the functions phi(t), psi(t), theta(t), f(t), etc,, then you need to provide additional equations.  In general, you need the number of equations to be equal to the numbers of unknowns.

@Fancypants In the worksheet you say that you expect partial_py1 to have terms involving the derivatives of S[y](t) with respect to py1.  But Maple hasn't been told that S[y](t) depends on py1, therefore it sets all those derivatives to zero.

There may be other errors in the worksheet but I don't quite see the aim of your calculations, and therefore I cannot suggest a remedy.

@spinoza Are you going to pass the solutions that you have obtained here as your own?  Or are you going to tell your teacher that Carl/Kitonum solved them for you?