@brian bovril Think i'll do that and see if i get the desired result. I have values for both  e and x. I don't know why i didn't think of fsolve initially. Thank you very much.

@Preben Alsholm  This looks good. I'll surely give it a trial. Thanks

@tomleslie  the result of convert(taylor(y(x+h), x=0, 3), polynom);

will not give me what i require. convert(taylor(y(x+h), h=0, 3), polynom) is what I want.

y(x) + y^'(x)*h + (y^''(x)*h²)/2! + (y^'''(x)*h³)/3!

Now, suppose i have a value for y(x), y^'''(x) and y^''(x) and h, i can i substitute it there...

@ThU It didn't work. if observe well, distributing the denominator to each expression in the numerator will bring out "cot" while the denominator cancels out on some terms.

@tomleslie I added the line Et:=Vector(N), where N:=5 at the start of the code but i noticed it wasn't working so I commented it out.

@tomleslie I noticed that even if I declared Et to be a Vector, Et[1] would still not work unless i make it double index like you've said. Any reason for that?

@Carl Love I understand. You've answered my next question. Thank you. I actually wrote a small piece of code to construct the matrix for me but later decided to ask the forum for alternative way(s) which is far far better.

@Rouben Rostamian  This is cool thanks it worked

@Carl Love My teacher gave a dazing and kinda ugly reply. Just joking please kindly explain

z:= Matrix(N-n,n):  Z:= Matrix(N,N-n):
<<A,z>|Z> + <Z|<z,B>>;

I have no sigle idea of what it interpretes

@Kitonum That takes care of the error. Thanks

@Kitonum I modified the code to look like this and still have the error

eul := proc(x,y)
    local a,b:
    a:=x:
    b:=y:
    for j from 1 to 3 do
          y[j]:= y + h*fxy(a,b):
          
          a:= 0 + j*h:
          b := y[j]:
          yt[j] := b:
    end do:
    return yt:
end proc:

@Carl Love I'm definitely dreaming in code. Nice one and thanks a million times.

@Carl Love not all the 11 points, connect the points (5 points to the left or 6 to the right )

(6*(20*z^5+137*z^4+450*z^3+850*z^2+900*z+420))/

(300*z^6-1490*z^5+4197*z^4-7800*z^3+9600*z^2-7200*z+2520).

I dont even mind if it is possible to shade the connected regions just like the sample graph i included.

@Carl Love Thanks alot. I'll try the code and try to take the modal value as advised.

@Carl Love This looks good. I'll definitely go and try it out then give you a reply. Thanks man.