Thanks for the suggestion, I have completed as per adviced.
My aviailable physical memory is 4GB 

Please have a look . I have asked some questions inside.
non_linear_P2_last_q.mw

I was trying some ways around and had some success but here is the sheet non_linear_P1.mw .

Please let me know in case of any ambiguity.

God bless you!!!

Is there any way to change the limits of integrals in the eq7?

I am interested in converting the limits to x[i]..x[i+1].

Though at this point, I am doing the following but I think there has to be some better way.

eq8:=eval(eq7, [a[j]= 0,  a[j+1/2]= 0,  a[j+1]= 0]):
eq9:=eq7-eq8=eq8;

eq10:=subs(phi[i](r)=phi[j](r), eq9):
collect(eq10,int);

eq11:=subs({int(diff(phi[j](r),r)*phi[j](r),r=1..delta)=

int( diff(phi[j](r),r)*phi[j](r),r=x[i]..x[i+1]),int(diff(phi[j](r),r)*phi[j+1](r),r=1..delta)

=int( diff(phi[j](r),r)*phi[j+1](r),r=x[i]..x[i+1]), int(diff(phi[j](r),r)*phi[j+1/2](r),

r=1..delta)=int( diff(phi[j](r),r)*phi[j+1/2](r),r=x[i]..x[i+1]),

 int(diff(phi[j+1/2](r),r)*phi[j](r),r=1..delta)=int( diff(phi[j+1/2](r),r)*phi[j](r),

r=x[i]..x[i+1]) }, eq10):

sorry therse x[i]'s has to be r[i]'s

It made me smile!!!

I think 'i' has to be unassigned, though I can see that implicitly definig 'i' would not be a problem, since,

my mesh in the latter case is ..., x[i]<x[i+1/2]<x[i+1], ...

,the integrals are for the ith element stiffness matrix, in the case of 1D problem with quadratic basis functions.

I have the following code which is working fine but I am not satisfied with the implementation, any suggestion to make it more aesthetic.

restart:
interface(rtablesize=infinity):
N:=3:
K:=Matrix(N):
x[i+1/2]:=x[i]+h/2:
x[i+1]:=x[i]+h:
phi[i]:=piecewise(t>=x[i] and t<=x[i+1],1-3*((t-x[i])/h)+2*(((t-x[i])^2)/h)):
phi[i+1/2]:=piecewise(t>=x[i] and t<=x[i+1], 1-4* (((t-x[i+1/2])^2)/h),0):
phi[i+1]:=piecewise(t>=x[i] and t<=x[i+1],1+3*(t-x[i+1])/h+2*((t-x[i+1])^2)/h,0):

K[1,1]:=int(diff(phi[i],t)*diff(phi[i],t), t=x[i]..x[i+1])assuming h>0:
K[2,2]:=simplify(int(diff(phi[i+1/2],t)*diff(phi[i+1/2],t), t=x[i]..x[i+1])) assuming h>0:
K[3,3]:=int(diff(phi[i+1],t)*diff(phi[i+1],t), t=x[i]..x[i+1]) assuming h>0:

K[1,2]:=simplify(int(diff(phi[i],t)*diff(phi[i+1/2],t), t=x[i]..x[i+1])) assuming h>0:
K[1,3]:=simplify(int(diff(phi[i],t)*diff(phi[i+1],t), t=x[i]..x[i+1])) assuming h>0:

K[2,1]:=simplify(int(diff(phi[i+1/2],t)*diff(phi[i],t), t=x[i]..x[i+1])) assuming h>0:
K[2,3]:=simplify(int(diff(phi[i+1/2],t)*diff(phi[i+1],t), t=x[i]..x[i+1])) assuming h>0:

K[3,1]:=simplify(int(diff(phi[i+1],t)*diff(phi[i],t), t=x[i]..x[i+1])) assuming h>0:
K[3,2]:=simplify(int(diff(phi[i+1],t)*diff(phi[i+1/2],t), t=x[i]..x[i+1])) assuming h>0:
K;

restart;
N:=10:
K:=Matrix(N-1):
for i from 1 to N do
 phi[i]:=piecewise(t>=x[i-1] and x[i]>t, (t-x[i])/(h), x[i]<t and t<=x[i+1], (x[i+1]-t)/(h), 0):
end do:
for i from 2 to N-2 do
    V1:=Int(diff(phi[i], t)*diff(phi[i], t), t=x[i]..x[i+1]):
    V2:=Int(diff(phi[i], t)*diff(phi[i+1], t), t=x[0]..x[N]):
    L1:=value(V1) assuming x[i-1]<x[i],x[i]<x[i+1];
    M1:=algsubs(x[i]-x[i-1]=h, L1);
    K[i,i]:=algsubs(x[i+1]-x[i]=h, M1);

    L2:=value(V2) assuming seq(x[i-1]<x[i], i=1..N);
    M2:=algsubs(x[i]-x[i-1]=h, L2);
    K[i,i+1]:=algsubs(x[i+1]-x[i]=h, M2);
    K[i-1,i]:=K[i,i+1];
 end do:
K;

So I have got this, its runing pretty as expected except the last statement

 > K[i-1,i]:=K[i,i+1].

There is just one problem that the algorithim is runing too slow even for this small value of N,

 since I am mainly interested in runing this for some large value of N, 

therefore I would really appreciate I someone could suggest some thing to make it fast.

The representation of phi[i]

phi[i]:=t->piecewise(t>=x[i-1] and x[i]>t, (t-x[i])/(h), x[i]<t and t<=x[i+1], (x[i+1]-t)/(h), 0);

produces the same output but in that case maple does not use heavside function but its definition.

phi[i]:=piecewise(t>=x[i-1] and x[i]>t, (t-x[i])/(h), x[i]<t and t<=x[i+1], (x[i+1]-t)/(h), 0);

@Carl Love 

Thanks for looking into it.

here is the command

phi[i]:=piecewise(t>=x[i-1] and x[i]>t, (t-x[i])/(h), x[i]

I tried it very extensively but you may agree that Maple online help browser isn't that good, I don't know whether they have done some thing in newer versions or not.

Afterwards after using bing.com, I came across some maple sheet by Dr. Robert Lopez(http://www.maplesoft.com/applications/view.aspx?SID=1742), the author is trying to compare ways to do by parts integration in maple, a nice overview of by parts integration in maple.

I have figure that out.

with(IntegrationTools):

integral_1:=  int( diff( u(x), x$2) * v(x), x=0..1):

Integral_2:= Parts(integral_1, v(x)) = D(u)(1)*v(1) - D(u)(0)*v(0) - int( diff(u(x), x) * diff(v(x), x), x=0..1)

subs({v(1) = 0, v(0)= 0}, Integral_2 ) = - int( diff(u(x), x) * diff(v(x), x), x=0..1)

This completes what I was tyring to accomplished.

I do have little point, I know define() is more powerful as compared to subs()  but I was unable to use that. Would some body help?

Thanks guys for kind words.

The two experts with complete different nature of expertise.

Thankyou very much!

Is it possible in maple to make the coordinate axes thick (or to rephrase my question, do we have a control on the thickness or colour of the coordinate axis )? as we do have a control on the thickness of the curve in maple.