Physics

% x'' + 2ζω x' + ω² x = 0, written as a first-order system y = [x; x']
w = 2*pi;
for zeta = [0.05 0.2 0.5 1]
  f = @(tt, y) [y(2); -2*zeta*w*y(2) - w^2*y(1)];
  [tt, y] = ode45(f, [0 4], [1; 0]);
  plot(tt, y(:, 1), 'LineWidth', 1.6, 'DisplayName', sprintf('ζ = %.2f', zeta))
  hold on
end
hold off

yline(0, 'k:')
legend show
title('Damped harmonic oscillator')
subtitle('ode45, adaptive Runge–Kutta')
xlabel('time (s)'), ylabel('displacement')

% u_t = α u_xx on [0, 1] with u = 0 at both ends, from a square pulse.
animate(8)
alpha = 0.02;
x = linspace(0, 1, 300);
u = zeros(size(x));
for k = 1:60                        % Fourier sine series solution
  bk = 2 * (cos(k*pi*0.35) - cos(k*pi*0.65)) / (k*pi);
  u = u + bk * sin(k*pi*x) * exp(-alpha * (k*pi)^2 * t);
end

area(x, u)
ylim([0 1.2])
title(sprintf('Heat diffusion   t = %.2f', t))
xlabel('position'), ylabel('temperature')

[TH, W] = meshgrid(linspace(-2*pi, 2*pi, 25), linspace(-3, 3, 19));
quiver(TH, W, W, -sin(TH))
hold on
for E = [-0.5 0.2 1 1.8]           % a few energy levels: ½ω² − cos θ = E
  th = linspace(-2*pi, 2*pi, 600);
  v = 2*(E + cos(th));  v(v < 0) = NaN;   % no motion where the energy is too low
  w = sqrt(v);
  plot(th, w, 'LineWidth', 1.5), plot(th, -w, 'LineWidth', 1.5)
end
hold off
xlim([-2*pi 2*pi]), ylim([-3 3])
title('Pendulum: θ'' = −sin θ'), xlabel('θ'), ylabel('ω')

[X, Y] = meshgrid(linspace(-3, 3, 160));
q = [1 -1]; px = [-1 1];
V = zeros(size(X));
for k = 1:2, V = V + q(k) ./ sqrt((X - px(k)).^2 + Y.^2 + 0.05); end
contourf(X, Y, max(min(V, 3), -3), 20)
colormap coolwarm
hold on
[Xq, Yq] = meshgrid(linspace(-3, 3, 17));
Ex = zeros(size(Xq)); Ey = Ex;
for k = 1:2
  d3 = ((Xq - px(k)).^2 + Yq.^2 + 0.05).^1.5;
  Ex = Ex + q(k) * (Xq - px(k)) ./ d3;
  Ey = Ey + q(k) * Yq ./ d3;
end
L = sqrt(Ex.^2 + Ey.^2);
quiver(Xq, Yq, Ex ./ L, Ey ./ L, 'Color', 'k')
hold off
axis equal
title('Dipole: potential and field direction')

animate(6)
[X, Y] = meshgrid(linspace(-6, 6, 220));
r1 = sqrt((X + 1.5).^2 + Y.^2);
r2 = sqrt((X - 1.5).^2 + Y.^2);
k = 2*pi/1.2; w = 2*pi/2;
u = cos(k*r1 - w*t) ./ sqrt(r1 + 1) + cos(k*r2 - w*t) ./ sqrt(r2 + 1);
imagesc(X(1, :), Y(:, 1), u)
axis xy, axis equal
colormap coolwarm
clim([-1.6 1.6])
title('Two point sources in phase')

x = linspace(0, 1, 400);
for n = 1:4
  E = n^2;                                        % energy in units of h²/(8mL²)
  psi = sqrt(2) * sin(n*pi*x);
  plot(x, E + 0.8*psi, 'LineWidth', 2), hold on
  plot(x, E + 0.4*psi.^2, ':')
  yline(E, 'k:')
end
hold off
xlabel('position x / L'), ylabel('energy (h²/8mL²)')
title('Wavefunctions (solid) and probability densities (dotted)')

g = 9.81; v0 = 50; ang = 45*pi/180;
k = 0.02;                                         % quadratic drag per unit mass (1/m)
f = @(tt, s) [s(3); s(4); -k*norm(s(3:4))*s(3); -g - k*norm(s(3:4))*s(4)];
[tt, s] = ode45(f, [0 7], [0; 0; v0*cos(ang); v0*sin(ang)]);
up = s(:, 2) >= 0;
plot(s(up, 1), s(up, 2), 'LineWidth', 2), hold on
T = 2*v0*sin(ang)/g; ts = linspace(0, T, 200);
plot(v0*cos(ang)*ts, v0*sin(ang)*ts - g*ts.^2/2, '--', 'LineWidth', 2)
hold off
axis equal
legend('with air drag', 'in vacuum')
xlabel('distance (m)'), ylabel('height (m)')
title(sprintf('Drag cuts the range from %.0f m to about %.0f m', v0^2*sin(2*ang)/g, max(s(up, 1))))

animate(4)
x = linspace(0, 1, 400);
w = 2*pi/4;                                       % fundamental: one cycle per loop
u = 0.8*sin(pi*x)*cos(w*t) + 0.35*sin(3*pi*x)*cos(3*w*t);
plot(x, u, 'LineWidth', 3), hold on
plot(x, 0.8*sin(pi*x) + 0.35*abs(sin(3*pi*x)), 'k:')
plot(x, -0.8*sin(pi*x) - 0.35*abs(sin(3*pi*x)), 'k:')
hold off
ylim([-1.3 1.3])
title('Modes 1 and 3 of a fixed string')
xlabel('position'), ylabel('displacement')