HOW DO I DO THAT IN MATLAB SERIES?

In this series, I am answering questions that students have asked me about MATLAB. Most of the questions relate to a mathematical procedure.

TOPIC
SUMMARY
INTRODUCTION
INPUTS
DISPLAYING INPUTS
THE CODE
DISPLAYING OUTPUTS

TOPIC

How do I solve a boundary value ordinary differential equation?

SUMMARY

% Language : Matlab 2008a;
% Authors : Autar Kaw;
% Mfile available at
% http://numericalmethods.eng.usf.edu/blog/ode_boundaryvalue.m;
% Last Revised : May 26, 2009;
% Abstract: This program shows you how to solve an
%           boundary value ordinary differential equation.
clc
clear all

INTRODUCTION

disp('ABSTRACT')
disp('   This program shows you how to solve')
disp('   a boundary value ordinary differential equation')
disp(' ')
disp('AUTHOR')
disp('   Autar K Kaw of http://autarkaw.wordpress.com')
disp(' ')
disp('MFILE SOURCE')
disp('   http://numericalmethods.eng.usf.edu/blog/ode_boundaryvalue.m')
disp(' ')
disp('LAST REVISED')
disp('   May 26, 2009')
disp(' ')

ABSTRACT
   This program shows you how to solve
   a boundary value ordinary differential equation
 
AUTHOR
   Autar K Kaw of http://autarkaw.wordpress.com
 
MFILE SOURCE
   http://numericalmethods.eng.usf.edu/blog/ode_boundaryvalue.m
 
LAST REVISED
   May 26, 2009

INPUTS

Solve the ordinary differential equation y"-0.0000002y=0.000000075*x*(75-x) Define x as a symbol

syms x
%The ODE
ode_eqn='D2y-0.000002*y=0.000000075*x*(75-x)';
% The initial conditions
iv_1='y(0)=0';
iv_2='y(75)=0';
% The value at which y is sought at
xval=37.5;

DISPLAYING INPUTS

disp('INPUTS')
func=['  The ODE to be solved is ' ode_eqn];
disp(func)
iv_explain=['  The boundary conditions are ' iv_1 '    ' iv_2];
disp(iv_explain)
fprintf('  The value of y is sought at x=%g',xval)
disp('  ')

INPUTS
  The ODE to be solved is D2y-0.000002*y=0.000000075*x*(75-x)
  The boundary conditions are y(0)=0    y(75)=0
  The value of y is sought at x=37.5

THE CODE

% Finding the solution of the ordinary differential equation
soln=dsolve(ode_eqn,iv_1,iv_2,'x');
% vpa below uses variable-precision arithmetic (VPA) to compute each
% element of soln to 5 decimal digits of accuracy
soln=vpa(soln,5);

DISPLAYING OUTPUTS

disp('  ')
disp('OUTPUTS')
output=['  The solution to the ODE is ' char(soln)];
disp(output)
value=subs(soln,x,xval);
fprintf('  The value of y at x=%g is %g',xval,value)
disp('  ')

  
OUTPUTS
  The solution to the ODE is -17756.*exp(.14142e-2*x)-19744.*exp(-.14142e-2*x)+37500.-2.8125*x+.37500e-1*x^2
  The value of y at x=37.5 is -0.0020632