Develop matlab function to solve problem

Complete and submit the following with your group. Write a user-defined MATLAB function that solves, with the shooting method in conjunction with the secant method, a second-order boundary value problem of the form: d^2y/dx^2 + f(x) dy/dx + g(x)y = h(x) for a lessthanorequalto x lessthanorequalto b with y(a) = Y_a and y(b) = Y_b For the function name and arguments use [x, y] = BVPShootSecant(fOFx, gOFx, hOFx, a, b, n, Ya, Yb, WL, WH). The input arguments fOFx, gOFx, and hOFx are names for the functions that calculate f(x), g(x), and h(x), respectively. They are dummy names for the anonymous or user-defined functions that are imported into BVPShootSecant. The arguments a and b define the domain of the solution, n is the number of sub-intervals, Ya and Yb are the boundary conditions, and WL and WH are the assumed slopes at x = a that are used in the first two solutions. Use ode45 to solve the ODEs. Iterate until the absolute value of the true error at x = b is smaller than 0.001. Test with problem #2 (parts a-b) of HW11, reproduced below: A spherical organism whose radius is R = 100 mu m consumes O_2 at a rate of k_O_2 = 25 s^-1 throughout its volume. The oxygen concentration at the surface of the organism is 3 times 10^22 m^-3, and its value in the interior of the organism is governed by the following one-dimensional diffusion equation: 1/r^2 d/dr (Dr^2 dn_O/dr) - k_O_1 n_O_2 = 0 where D = 10^-8 m^2/s is the diffusivity of O_2. By symmetry, the boundary condition at r = 0 is dn_O/dr|_r = 0 = 0. (a) Rewrite the ODE as a system of two first-order ODEs and write down the boundary conditions. (b) Solve the system of first-order ODEs using the MATLAB built-in function bvp4c for the steady state n_O_2(r). Note that r = 0 appears as a singularity on the right-hand side of one of the equations in part (a), and must therefore be replaced with a small (compared to R), non-zero value, e g. 10^-2 m. Solve the equations from r = 10^-7 m to r = R by dividing the domain into 100 subintervals, and using initial guesses of n_O_2 = 3 times 10^22 and dn_O_2/dr = 0. Make a plot of n_O_2 versus r