Matlab may be used to solve differential equations symbolically as well as numerically. Well, actually Maple is being used behind the scenes. This can be seen whenever an error occurs with the solver 'Error using ==> maple at 129'. Regardless, say you wanted to solve for a spring mass damper system under a forced oscillation:
m*d(d(x))+c*d(x)+k*(x)=f*cos(w*t)
The solution is lots of fun to do by hand, but faster and easier in Matlab using dsolve. Dsolve calls Maple to symbolically solve the system. The solver can solve for a single equation with boundary conditions or for systems of equations. For this case its one second order non harmonic with two initial conditions (I will set the initial position to 0.05 to keep the solution 'interesting').
dsolve('m*D2x+c*Dx+k*x=f*cos(w*t)','x(0)=0.05','Dx(0)=0')
Running this returns:
1/40*exp(-1/2*(c-(c^2-4*k*m)^(1/2))/m*t)*(2*w^2*m^2*k-2*k^2*m+40*m*f*k+k*c^2-(c^2-4*k*m)^(1/2)*c*k+20*(c^2-4*k*m)^(1/2)*c*f-20*f*c^2)*(-c*(c^2-4*k*m)^(1/2)-c^2+4*k*m)/(-2*w^2*m^2*c^2+8*w^2*m^3*k-c^4-8*k^2*m^2+6*c^2*k*m+c^3*(c^2-4*k*m)^(1/2)-4*m*(c^2-4*k*m)^(1/2)*c*k)/k+1/40*exp(-1/2*(c+(c^2-4*k*m)^(1/2))/m*t)*(2*w^2*m^2*k-2*k^2*m+40*m*f*k+k*c^2+(c^2-4*k*m)^(1/2)*c*k-20*(c^2-4*k*m)^(1/2)*c*f-20*f*c^2)*(4*k*m-c^2+c*(c^2-4*k*m)^(1/2))/(2*m^2*w^2+c^2-2*k*m+c*(c^2-4*k*m)^(1/2))/k/(-c^2+4*k*m)-f*(-cos(w*t)*k+cos(w*t)*w^2*m-sin(w*t)*w*c)/(c^2*w^2+k^2-2*k*m*w^2+w^4*m^2)
Fun, but say I have values for m, c, k, f and w…
subs(subs(subs(subs(subs(dsolve('m*D2x+c*Dx+k*x=f*cos(w*t)','x(0)=0.05','Dx(0)=0'),'k',2000),'c',20'),'m',5),'w',pi),'f',10)
Running this returns:
1/80000*exp(-1/10*(20-(-39600)^(1/2))*t)*(100000*pi^2-35280000-36000*(-39600)^(1/2))*(39600-20*(-39600)^(1/2))/(1980000*pi^2-776160000-792000*(-39600)^(1/2))+1/3168000000*exp(-1/10*(20+(-39600)^(1/2))*t)*(100000*pi^2-35280000+36000*(-39600)^(1/2))*(20*(-39600)^(1/2)+39600)/(50*pi^2-19600+20*(-39600)^(1/2))-10*(-2000*cos(pi*t)+5*cos(pi*t)*pi^2-20*sin(pi*t)*pi)/(-19600*pi^2+4000000+25*pi^4)
This is better, but it could be simpler…
vpa(subs(subs(subs(subs(subs(dsolve('m*D2x+c*Dx+k*x=f*cos(w*t)','x(0)=0.05','Dx(0)=0'),'k',2000),'c',20'),'m',5),'w',pi),'f',10))
Running this after setting the digits to five (by running 'digits 5') returns:
(.22439e-1-.22422e-2*i)*exp((-2.0000+19.900*i)*t)+(.22439e-1+.22420e-2*i)*exp((-2.0000-19.900*i)*t)+.51214e-2*cos(3.1416*t)+.16496e-3*sin(3.1416*t)
This is much more manageable. I wonder what it looks like:
plot([0:.01:1], subs(vpa(subs(subs(subs(subs(subs(dsolve('m*D2x+c*Dx+k*x=f*cos(w*t)','x(0)=0.05','Dx(0)=0'),'k',2000),'c',20'),'m',5),'w',pi),'f',10)),'t',[0:.01:1]))
Running this plots:
Note that there are better ways to plot this. For example, one could take the simplified expression and make it into an expression and use fplot. This will plot faster than substituting for each data point and evaluating the string.