% York U PHYS 2030 (09.15.14) % use basic Euler method to solve the following differential equation % f'(x)= x^2 + 1 clear % ************************ % User Inputs Fprime = @(x)(x.^2+1); % function describing the ODE Fsol= @(x,IC)(1/3*x.^3 + x + IC); % solution (known; IC is the initial condition, f0) stepsize= 0.5; f0= 0.0; % initial condition at xI [i.e. f(xI)] xB= [0 2]; % interval to solve over (i.e., boundaries) % ************************ F(1)= f0; %initialize first value k=1; % counter for j=xB(1):stepsize:xB(2) if (j == xB(1)), F(k)= f0; % handle initial condition else F(k)= F(k-1) + stepsize*(Fprime(j)); % apply Euler's method end x(k)= j; % keep track of x for plotting k= k+ 1; % increment counter end % visualize figure(1); clf; plot(x,F, 'o--','LineWidth',2); grid on; hold on xlabel('x'); ylabel('f(x)'); title('Solution to df/dx= x^2 +1') % now plot exact solution and derivative plot(x,Fsol(x,f0), 'r.-','LineWidth',1); plot(x,Fprime(x), 'kd-'); legend('f(x) (numeric solution)', 'f(x) (true solution)', 'df/dx (from ODE)','Location','NorthWest')