I am trying to solve a system of 4 odes by using the following program:

with(plots):

Pr:=6.2:Ha:=10^(-12):Nt:=0.1:Nb:=0.1:Nr:=0.1:Gr:=0.1:Ec:=0.1:Rb:=0.1:Lb:=1:Pe:=1:tw:=0.2:Le:=5:

fcns:={f(y),T(y),S(y),G(y)}:

sys:=diff(f(y),y$3)+f(y)*diff(f(y),y$2)-diff(f(y),y)*diff(f(y),y)-Ha*diff(f(y),y)+Gr*(T(y)-Nr*S(y)-Rb*G(y))=0,

diff(T(y),y$2)+Pr*f(y)*diff(T(y),y)+Nb*diff(S(y),y)*diff(T(y),y)+Nt*diff(T(y),y$2)*diff(T(y),y$2)+Pr*Ec*(diff(f(y...