Question 24.8: Use bvp4c to solve the following second-order ODE d²y/dx² + ...
Use bvp4c to solve the following second-order ODE
\frac{d^2 y}{d x^ 2} + y = 1
subject to the boundary conditions
y(0) = 1
\\ y(π ∕ 2) = 0
First, express the second-order equation as a pair of first-order ODEs
\frac{d y}{dx} = z
\frac{dz}{dx} = 1 − y
Next, set up a function to hold the first-order ODEs
function dy = odes(x,y)
dy = [y(2); 1 − y(1)];
We can now develop the function to hold the boundary conditions. This is done just like a roots problem in that we set up two functions that should be zero when the boundary conditions are satisfied. To do this, the vectors of unknowns at the left and right boundaries are defined as ya and yb. Hence, the first condition, y(0) = 1, can be formulated as ya(1) − 1; whereas the second condition, y(π/2) = 0, corresponds to yb(1).
function r = bcs(ya,yb)
r = [ya(1) − 1; yb(1)];
Finally, we can set up solinit to hold the initial mesh and solution guesses with the bvpinit function. We will arbitrarily select 10 equally spaced mesh points, and initial guesses of y = 1 and z = dy/dx = –1.
solinit = bvpinit (linspace(0,pi/2,10),[1,−1]);
The entire script to generate the solution is
clc
solinit = bvpinit(linspace(0,pi/2,10),[1,−1]);
sol = bvp4c(@odes,@bcs,solinit);
x = linspace(0,pi/2);
y = deval(sol,x);
plot(x,y(1,:))
where deval is a built-in MATLAB function which evaluates the solution of a differential equation problem with the general syntax
yxint = deval(sol,xint)
where deval evaluates the solution at all the values of the vector xint, and sol is the structure returned by the ODE problem solver (in this case, bvp4c).
When the script is run the plot below is generated. Note that the script and functions developed in this example can be applied to other boundary value problems with minor modifications. Several end-of-chapter problems are included to test your ability to do just that.
