> 
# # Define some parameters # sigma := 10; N := 0; beta := 1; alpha := 1; PDE1 := diff(w(X, theta, t), X, X, X, X)+2*alpha^2*(diff(w(X, theta, t), theta, theta, X, X))+alpha^4*(diff(w(X, theta, t), theta, theta, theta, theta))N*(diff(w(X, theta, t), X, X))+diff(w(X, theta, t), t, t)beta*w(X, theta, t)sigma = 0


(1) 
> 
# # Define the PDES # PDEs:= { diff(w(X, theta, t), X, X, X, X)+2*alpha^2*(diff(w(X, theta, t), theta, theta, X, X))+alpha^4*(diff(w(X, theta, t), theta, theta, theta, theta))N*(diff(w(X, theta, t), X, X))+diff(w(X, theta, t), t, t)beta*w(X, theta, t)sigma = 0 };


(2) 
> 
# # Set of boundary conditions at x=1. # bcs1:= { D[1](w)(1,theta, t) = 0, w(1,theta, t) = 0 };


(3) 
> 
# # Set of boundary conditions at x=0 # bcs2:= { w(0,theta, t)=0, D[1](w)(0,theta, t)=0 };


(4) 
> 
# # Set of boundary conditions at t=0 # bcs3:= { w(x,theta,0)=0, D[2](w)(x,theta,0)=0 };


(5) 
> 
pdsolve( PDEs, `union`(bcs1, bcs2, bcs3), numeric);

