Calculus

% ─── Numerical exploration of lim sin(x)/x as x→0 ───────────────
x_vals = [0.5, 0.1, 0.01, 1e-4, 1e-8];
printf("x           sin(x)/x\n");
for x = x_vals
  printf("%.2e    %.15f\n", x, sin(x)/x);
endfor
% Catastrophic cancellation demo — do NOT use tiny h for this limit
printf("\nh = 1e-16: sin(h)/h = %.15f  (numerical noise!)\n", sin(1e-16)/1e-16);

% ─── ε-δ verification: lim_{x→3}(2x−1) = 5 ─────────────────────
a = 3;  L = 5;  epsilon = 0.01;  delta = epsilon/2;
xs = a + delta*[-0.99, -0.5, 0.5, 0.99];
printf("\nε=%.4f δ=%.4f verification:\n", epsilon, delta);
for x = xs
  fx = 2*x - 1;
  printf("  x=%.4f  f=%.4f  |f-L|=%.4f  within_ε=%d\n", ...
         x, fx, abs(fx-L), abs(fx-L) < epsilon);
endfor

% ─── L'Hôpital: lim sin(x)/x via differentiation ─────────────────
% Numerator limit: d/dx sin(x) at 0 = cos(0) = 1
% Denominator lim: d/dx x at 0 = 1  →  ratio = 1
printf("\nL'Hopital result: cos(0)/1 = %.15f\n", cos(0));

% ─── IVT: bisection root of f(x)=x³−x−2 on [1,2] ────────────────
function c = bisect(f, a, b, tol)
  while (b - a)/2 > tol
    m = (a + b)/2;
    if f(m) == 0; c = m; return; endif
    if sign(f(m)) == sign(f(a)); a = m; else; b = m; endif
  endwhile
  c = (a + b)/2;
endfunction

f    = @(x) x^3 - x - 2;
root = bisect(f, 1, 2, 1e-12);
printf("\nBisection root of x³−x−2 = %.14f  f(root)=%.2e\n", root, f(root));

% ─── Continuity check: is |x|/x continuous at 0? ─────────────────
for x = [-0.01, -1e-10, 0, 1e-10, 0.01]
  if x == 0; val = 'undef'; else; val = num2str(x/abs(x)); endif
  printf("  sign(%.2e) = %s\n", x, val);
endfor
printf("→ Limits from left/right disagree: NOT continuous at 0.\n");

% ─── Complex-step derivative (machine-precision, no cancellation) ─
%   f′(x) ≈ Im[f(x+ih)]/h  for tiny h — provably exact to h²
function d = csd(f, x)
  h = 1e-200;
  d = imag(f(x + 1i*h)) / h;
endfunction

printf("Complex-step derivative verification:\n");
tests = {'xⁿ',  @(x) x^5,         @(x) 5*x^4;
         'eˣ',   @exp,            @exp;
         'sin x',@sin,            @cos;
         'ln x', @log,            @(x)1/x;
         'xln x',@(x)x*log(x),   @(x)log(x)+1};
x0 = 1.7;
for i = 1:rows(tests)
  ex = tests{i,3}(x0);
  nu = csd(tests{i,2}, x0);
  printf("  d/dx %-8s exact=%+.10f  csd=%+.10f  err=%.1e\n", ...
         tests{i,1}, ex, nu, abs(ex-nu));
endfor

% ─── Central difference vs complex step accuracy comparison ───────
f  = @(x) sin(x.^2) .* exp(-x);
fp = @(x) (2*x.*cos(x^2) - sin(x^2)) .* exp(-x);
x0 = pi/3;  exact = fp(x0);
printf("\nError vs step size h:\n");
for h = [1e-1,1e-3,1e-6,1e-9,1e-12]
  cd_err = abs((f(x0+h)-f(x0-h))/(2*h) - exact);
  cs_err = abs(csd(@(x)f(x), x0) - exact);  % always ~1e-16 regardless of h
  printf("  h=%.0e  central_diff err=%.1e\n", h, cd_err);
endfor
printf("  complex-step err = %.1e (independent of h)\n", csd(@(x)f(x),x0)-exact);

% ─── Chain rule: d/dx sin(x²) = 2x·cos(x²) ──────────────────────
g  = @(x) sin(x.^2);
gp = @(x) 2*x.*cos(x.^2);
for x = [0.5, 1, pi/4]
  printf("  d/dx sin(x²)|x=%.3f: exact=%+.10f  csd=%+.10f\n", x, gp(x), csd(g,x));
endfor

% ─── Implicit diff: slope of ellipse 3x²+2y²=5 at (1,1) ──────────
% dy/dx = −3x/(2y)
x0=1; y0=1;
slope = -3*x0/(2*y0);
printf("\nImplicit: slope of 3x²+2y²=5 at (1,1): dy/dx = %.6f\n", slope);

% ─── Verify all 6 trig derivatives at x = π/5 ────────────────────
function d = csd(f, x); d = imag(f(x+1e-200i))*1e200; endfunction

x = pi/5;
rules = {'sin x',  @sin,  @cos;
         'cos x',  @cos,  @(x)-sin(x);
         'tan x',  @tan,  @(x)1/cos(x)^2;
         'csc x',  @csc,  @(x)-csc(x)*cot(x);
         'sec x',  @sec,  @(x)sec(x)*tan(x);
         'cot x',  @cot,  @(x)-csc(x)^2};
printf("Trig derivative verification at x=π/5:\n");
for i=1:rows(rules)
  ex = rules{i,3}(x);
  nu = csd(rules{i,2}, x);
  printf("  d/dx %-7s exact=%+.10f  csd=%+.10f  err=%.1e\n", ...
         rules{i,1}, ex, nu, abs(ex-nu));
endfor

% ─── Verify inverse trig derivatives ──────────────────────────────
printf("\nInverse trig derivatives:\n");
iv = {'arcsin', @asin, @(x)1/sqrt(1-x^2);
      'arccos', @acos, @(x)-1/sqrt(1-x^2);
      'arctan', @atan, @(x)1/(1+x^2)};
x2 = 0.4;
for i=1:3
  ex = iv{i,3}(x2);
  nu = csd(iv{i,2}, x2);
  printf("  d/dx %-8s at 0.4: exact=%+.10f  csd=%+.10f  err=%.1e\n", ...
         iv{i,1}, ex, nu, abs(ex-nu));
endfor

% ─── nth derivative of sin: d^n/dx^n sin(x) = sin(x + n*π/2) ─────
x0 = pi/3;
printf("\nnth derivative of sin x at x=π/3:\n");
for n = 0:7
  val = sin(x0 + n*pi/2);
  printf("  n=%d: sin(x+nπ/2) = %+.8f\n", n, val);
endfor

% ─── Tricky chain rule: d/dx sin(cos(x)) ─────────────────────────
f  = @(x)sin(cos(x));
fp = @(x)cos(cos(x)).*(-sin(x));
x_grid = linspace(0,2*pi,5);
printf("\nd/dx sin(cos(x)):\n");
for xv = x_grid
  printf("  x=%.4f  exact=%+.8f  csd=%+.8f\n", xv, fp(xv), csd(f,xv));
endfor

% ─── Plot: sin, cos and their derivatives ─────────────────────────
x = linspace(-2*pi, 2*pi, 1000);
figure('Name','Trig Derivatives');
plot(x,sin(x),'b-','lw',2, x,cos(x),'r--','lw',2);
legend('sin x',"d/dx sin x = cos x");
grid on; title('sin x and its derivative');

% ─── Newton-Raphson with convergence tracking ─────────────────────
function [root,iters,errs] = newton(f, fp, x0, tol)
  x = x0; errs = [];
  for k = 1:50
    fx = f(x);  errs(end+1) = abs(fx);
    if abs(fx) < tol; root=x; iters=k; return; endif
    x = x - fx/fp(x);
  endfor
  root=x; iters=50;
endfunction

% Find √2: solve x²−2=0
f = @(x)x^2-2;  fp = @(x)2*x;
[r,k,e] = newton(f, fp, 1.5, 1e-15);
printf("√2 by Newton-Raphson:\n");
for i=1:length(e)
  printf("  iter %d: |f|=%.3e  (digits ≈ %.1f)\n", i, e(i), -log10(e(i)+1e-16));
endfor
printf("  Result: %.15f  error=%.2e\n", r, abs(r-sqrt(2)));

% Quadratic convergence: each step squares the error
printf("\nConvergence ratio |eₙ₊₁|/|eₙ|²:\n");
for i=1:length(e)-2
  printf("  step %d→%d: %.4f (should→constant)\n", i, i+1, e(i+1)/e(i)^2);
endfor

% ─── Optimize f(x) = x⁴−4x²+x using Newton on f′=0 ──────────────
f2  = @(x) x^4 - 4*x^2 + x;
fp2 = @(x) 4*x^3 - 8*x + 1;
fpp = @(x) 12*x^2 - 8;
printf("\nCritical points of x⁴−4x²+x:\n");
for x0 = [-1.5, 0.1, 1.6]
  [r,~,~] = newton(fp2, fpp, x0, 1e-13);
  typ = '?'; if fpp(r)>0; typ='min'; elseif fpp(r)<0; typ='max'; endif
  printf("  x₀=%.2f → x*=%.8f  f=%.6f  f''=%.4f → %s\n", x0,r,f2(r),fpp(r),typ);
endfor

% ─── Linearisation: approximate sin(0.1) near a=0 ────────────────
a = 0;
printf("\nLinearisation of sin(x) near x=0:\n");
for h = [0.1, 0.3, 0.5, 1.0]
  lin   = sin(a) + cos(a)*(h - a);   % = h (since sin(0)=0, cos(0)=1)
  exact = sin(h);
  printf("  x=%.1f: linear=%.6f  exact=%.6f  err=%.2e\n", h, lin, exact, abs(lin-exact));
endfor

% ─── Riemann sums: left, midpoint, right ─────────────────────────
function s = riemann(f, a, b, n, rule)
  h  = (b-a)/n;
  switch rule
    case 'left';  xs = a + (0:n-1)*h;
    case 'right'; xs = a + (1:n)*h;
    case 'mid';   xs = a + (0.5:1:n-0.5)*h;
  endswitch
  s = h * sum(f(xs));
endfunction

f     = @(x) sin(x) + 1.2;
a=0; b=pi;
exact = quadgk(f, a, b);
printf("∫₀^π (sin x + 1.2) dx = %.15f\n", exact);
printf("\nRiemann sum convergence:\n");
for n = [5, 20, 100, 1000]
  L = riemann(f,a,b,n,'left');
  M = riemann(f,a,b,n,'mid');
  printf("  n=%4d: left=%.8f  mid=%.8f  mid_err=%.2e\n", n,L,M,abs(M-exact));
endfor

% ─── Composite trapezoid and Simpson's ───────────────────────────
function s = trapz_rule(f,a,b,n)
  x = linspace(a,b,n+1); h=(b-a)/n;
  s = h/2*(f(x(1))+2*sum(f(x(2:end-1)))+f(x(end)));
endfunction
function s = simpsons(f,a,b,n)  % n must be even
  x = linspace(a,b,n+1); h=(b-a)/n;
  s = h/3*(f(x(1)) + 4*sum(f(x(2:2:end-1))) + 2*sum(f(x(3:2:end-2))) + f(x(end)));
endfunction

g = @(x) exp(-x.^2/2) / sqrt(2*pi);   % Gaussian PDF
exact_g = quadgk(g, -3, 3);
printf("\nGaussian CDF P(|X|<3) via quadrature:\n");
printf("  exact = %.15f\n", exact_g);
printf("  trapz (n=100) err = %.2e\n", abs(trapz_rule(g,-3,3,100)-exact_g));
printf("  Simp  (n=100) err = %.2e\n", abs(simpsons(g,-3,3,100)-exact_g));

% ─── FTC Part 1 demo: d/dx ∫₀ˣ sin(t²) dt ≈ sin(x²) ─────────────
F  = @(x) quadgk(@(t)sin(t.^2), 0, x);
Fp = @(x) sin(x^2);   % FTC Part 1 prediction
h  = 1e-5;
x0 = 1.3;
numerical_deriv = (F(x0+h) - F(x0-h))/(2*h);
printf("\nFTC Part 1: d/dx ∫₀ˣ sin(t²)dt at x=1.3:\n");
printf("  FTC prediction sin(x²) = %.10f\n", Fp(x0));
printf("  Numerical deriv of F(x) = %.10f\n", numerical_deriv);

% ─── Verify all basic trig integrals ─────────────────────────────
tests = {'∫ sin',     @sin,           @(x)-cos(x);
         '∫ cos',     @cos,           @sin;
         '∫ sec²',    @(x)sec(x).^2,  @tan;
         '∫ sec·tan', @(x)sec(x).*tan(x), @sec};
a=0.3; b=1.1;
printf("Trig integral verification [%.1f, %.1f]:\n", a, b);
for i=1:rows(tests)
  ftc  = tests{i,3}(b) - tests{i,3}(a);
  num  = quadgk(tests{i,2}, a, b);
  printf("  %-12s FTC=%.10f  quad=%.10f  err=%.1e\n", ...
         tests{i,1}, ftc, num, abs(ftc-num));
endfor

% ─── Fourier orthogonality: ∫₀^2π sin(mx)sin(nx) dx ─────────────
printf("\nFourier inner products (should be π on diagonal, 0 off):\n");
for m=1:4
  for n=1:4
    v = quadgk(@(x)sin(m*x).*sin(n*x), 0, 2*pi);
    printf("%6.3f", v/pi);
  endfor; printf("\n");
endfor

% ─── Gaussian integral: ∫_{-∞}^{∞} e^{-x²} dx = √π ─────────────
I = quadgk(@(x)exp(-x.^2), -30, 30);
printf("\n∫e^(-x²)dx = %.15f  √π = %.15f\n", I, sqrt(pi));

% ─── Power reduction: ∫₀^π sin⁴(x)dx = 3π/8 ────────────────────
I4 = quadgk(@(x)sin(x).^4, 0, pi);
printf("∫₀^π sin⁴x dx = %.12f  (3π/8=%.12f)\n", I4, 3*pi/8);

% ─── Trig substitution: ∫₀^1 √(1−x²)dx = π/4 (unit semicircle) ─
Iq = quadgk(@(x)sqrt(1-x.^2), 0, 1);
printf("∫₀^1 √(1−x²)dx = %.15f  (π/4=%.15f)\n", Iq, pi/4);

% ─── Wallis product: ∫₀^{π/2} sinⁿ dx via reduction formula ──────
function W = wallis(n)
  if n==0; W=pi/2; elseif n==1; W=1; else; W=(n-1)/n*wallis(n-2); endif
endfunction
printf("\nWallis integrals ∫₀^{π/2} sinⁿ(x) dx:\n");
for n=0:6
  printf("  n=%d: formula=%.8f  quad=%.8f\n", n, wallis(n), ...
         quadgk(@(x)sin(x).^n, 0, pi/2));
endfor

% ─── IBP: ∫ x·sin(x) dx = −x·cos(x) + sin(x) ────────────────────
antideriv = @(x) -x.*cos(x) + sin(x);
printf("∫₀^π x·sin(x) dx:\n");
printf("  FTC: %.15f  (exact: π=%.15f)\n", antideriv(pi)-antideriv(0), pi);
printf("  quad: %.15f\n", quadgk(@(x)x.*sin(x), 0, pi));

% ─── IBP: ∫ x²·eˣ dx (needs IBP twice) ──────────────────────────
% Antideriv: eˣ(x²−2x+2)
A2 = @(x) exp(x).*(x.^2 - 2*x + 2);
printf("\n∫₀^1 x²·eˣ dx:\n");
printf("  FTC: %.15f\n", A2(1)-A2(0));
printf("  quad: %.15f\n", quadgk(@(x)x.^2.*exp(x), 0, 1));

% ─── Reduction formula verification ──────────────────────────────
function W = wallis(n)
  if n==0; W=pi/2; elseif n==1; W=1; else; W=(n-1)/n*wallis(n-2); endif
endfunction
printf("\nReduction formula ∫₀^{π/2} sinⁿx dx:\n");
for n = [2,4,6,8]
  w = wallis(n);
  q = quadgk(@(x)sin(x).^n, 0, pi/2);
  printf("  n=%d: formula=%.10f  quad=%.10f  err=%.1e\n", n, w, q, abs(w-q));
endfor

% ─── Partial fractions: ∫ 1/(x²−1) dx ───────────────────────────
% 1/(x²−1) = (1/2)/(x−1) − (1/2)/(x+1)
% Antideriv = (1/2)·ln|(x−1)/(x+1)|
A3 = @(x) 0.5*log(abs((x-1)./(x+1)));
I_pf  = quadgk(@(x)1./(x.^2-1), 2, 5);
I_ftc = A3(5) - A3(2);
printf("\n∫₂^5 1/(x²−1) dx: quad=%.12f  FTC=%.12f  err=%.1e\n", I_pf, I_ftc, abs(I_pf-I_ftc));

% ─── u-sub: ∫ x·√(x²+1) dx = (x²+1)^{3/2}/3 ────────────────────
A4 = @(x) (x.^2+1).^(3/2)/3;
printf("∫₀^2 x√(x²+1) dx: FTC=%.10f  quad=%.10f\n", ...
       A4(2)-A4(0), quadgk(@(x)x.*sqrt(x.^2+1), 0, 2));

% ─── Taylor series for sin x: convergence ────────────────────────
function s = taylor_sin(x, N)
  s = 0;
  for n = 0:N-1
    s = s + (-1)^n * x^(2*n+1) / factorial(2*n+1);
  endfor
endfunction

x = pi/3;
printf("Taylor sin(π/3) convergence:\n");
for N = 1:9
  ap = taylor_sin(x, N);
  printf("  N=%d: %.15f  err=%.2e\n", N, ap, abs(ap-sin(x)));
endfor

% ─── Euler's formula: e^{iθ} = cos θ + i sin θ ───────────────────
theta = pi/5;
lhs = exp(1i*theta);
rhs = cos(theta) + 1i*sin(theta);
printf("\nEuler e^{iπ/5}:  LHS=%.10f%+.10fi\n", real(lhs), imag(lhs));
printf("               RHS=%.10f%+.10fi\n", real(rhs), imag(rhs));
printf("  e^{iπ}+1 = %.2e + %.2ei (should be 0)\n", ...
       real(exp(1i*pi)+1), imag(exp(1i*pi)+1));

% ─── Leibniz series and faster Machin formula for π ──────────────
printf("\nLeibniz π/4 partial sums:\n");
s = 0;
for n=0:9
  s = s + (-1)^n/(2*n+1);
  printf("  n=%d: 4s=%.10f  err=%.2e\n", n, 4*s, abs(4*s-pi));
endfor
% Machin: π/4 = 4·atan(1/5) − atan(1/239) — converges far faster
pi_machin = 4*(4*atan(1/5) - atan(1/239));
printf("Machin π=%.15f  err=%.2e\n", pi_machin, abs(pi_machin-pi));

% ─── Taylor coefficients via FFT (Cauchy's coefficient formula) ──
function c = taylor_fft(f, a, n, r)
  % c(k) = f^(k)(a)/k! via |f(a+r·e^{2πij/N})|/N
  N  = max(256, 2^ceil(log2(n+1)));
  th = 2*pi*(0:N-1)/N;
  z  = a + r*exp(1i*th);
  fv = arrayfun(f, z);
  cf = real(ifft(fv));
  c  = cf(1:n+1) ./ (r.^(0:n));
endfunction

c = taylor_fft(@sin, 0, 8, 1);
printf("\nTaylor coefficients of sin(x) at 0 (via FFT):\n");
for k=0:8
  if abs(c(k+1))>1e-10
    printf("  k=%d: %.10f  (1/k!=%.10f)\n", k, c(k+1), (-1)^((k-1)/2)/factorial(k));
  endif
endfor

% ─── Exponential decay: y'=−0.5y, y(0)=3 ────────────────────────
[t,y] = ode45(@(t,y)-0.5*y, [0 10], 3);
exact = @(t) 3*exp(-0.5*t);
err   = max(abs(y - exact(t)));
printf("Decay y'=−0.5y: max ode45 err = %.2e\n", err);

% ─── SHO: y''+4y=0, y(0)=1, y'(0)=0 → y=cos(2t) ─────────────────
sho = @(t,y)[y(2); -4*y(1)];
[ts,ys] = ode45(sho, [0 4*pi], [1;0]);
err_sho = max(abs(ys(:,1) - cos(2*ts)));
printf("SHO y''+4y=0: max err vs cos(2t) = %.2e\n", err_sho);

% ─── Damped oscillator: y''+2cy'+4y=0 for c=0.3,2,4 ─────────────
printf("\nDamped oscillator at t=10:\n");
for c = [0.3, 2, 4]
  ode = @(t,y)[y(2); -2*c*y(2)-4*y(1)];
  [td,yd] = ode45(ode,[0 10],[1;0]);
  disc = (2*c)^2 - 16;
  if disc<0; typ='under'; elseif disc==0; typ='crit'; else; typ='over'; endif
  printf("  c=%.1f (%sdamped): y(10)=%.6f\n", c, typ, interp1(td,yd(:,1),10));
endfor

% ─── Logistic equation ────────────────────────────────────────────
r=1.5; K=100; P0=5;
logistic_exact = @(t) K ./ (1 + (K-P0)/P0 * exp(-r*t));
[tl,Pl] = ode45(@(t,P)r*P*(1-P/K), [0 8], P0);
err_log = max(abs(Pl - logistic_exact(tl)));
printf("\nLogistic: P(8)=%.4f  exact=%.4f  err=%.2e\n", ...
       Pl(end), logistic_exact(8), err_log);

% ─── Lotka-Volterra (predator-prey) ──────────────────────────────
a=1.5; b=1; d=3; g=1;
lv = @(t,y)[a*y(1)-b*y(1)*y(2); d*y(1)*y(2)-g*y(2)];
[tlv,ylv] = ode45(lv,[0 15],[1;0.5]);
figure('Name','Predator-Prey');
plot(tlv,ylv(:,1),'b-','lw',2, tlv,ylv(:,2),'r-','lw',2);
legend('Prey','Predator'); grid on; title('Lotka-Volterra equations');

% ─── Partial derivatives via complex step ────────────────────────
f  = @(x,y) sin(x.*y) + x.^2.*cos(y);
h  = 1e-200;
x0 = pi/4;  y0 = pi/3;
dfdx_num = imag(f(x0+1i*h, y0))/h;
dfdy_num = imag(f(x0, y0+1i*h))/h;
% Exact: ∂f/∂x = y·cos(xy)+2x·cos(y),  ∂f/∂y = x·cos(xy)−x²·sin(y)
dfdx_ex = y0*cos(x0*y0)+2*x0*cos(y0);
dfdy_ex = x0*cos(x0*y0)-x0^2*sin(y0);
printf("∂f/∂x: exact=%.12f  cs=%.12f  err=%.1e\n", dfdx_ex, dfdx_num, abs(dfdx_ex-dfdx_num));
printf("∂f/∂y: exact=%.12f  cs=%.12f  err=%.1e\n", dfdy_ex, dfdy_num, abs(dfdy_ex-dfdy_num));

% ─── Gradient descent on f(x,y) = x²+2y²−xy+x ───────────────────
f2   = @(v) v(1)^2 + 2*v(2)^2 - v(1)*v(2) + v(1);
grad = @(v) [2*v(1)-v(2)+1; 4*v(2)-v(1)];
v = [4;3];  lr = 0.08;
printf("\nGradient descent (lr=0.08):\n");
for k=1:40
  v = v - lr*grad(v);
  if mod(k,8)==0
    printf("  k=%2d: v=[%.6f,%.6f]  f=%.8f  |∇|=%.2e\n",...
           k, v(1),v(2), f2(v), norm(grad(v)));
  endif
endfor
% Exact: ∇f=0 → 2x−y+1=0, 4y−x=0 → x=−4/7, y=−1/7
printf("  Exact min: x=%.8f, y=%.8f\n", -4/7, -1/7);

% ─── Hessian test for f(x,y)=x³−3x+y³−3y ────────────────────────
% Critical pts: 3x²−3=0 → x=±1; 3y²−3=0 → y=±1
printf("\nHessian classification of x³−3x+y³−3y:\n");
for cx = [-1,1]
  for cy = [-1,1]
    fxx=6*cx; fyy=6*cy; fxy=0;
    D = fxx*fyy - fxy^2;
    if     D>0&&fxx>0; t='min';
    elseif D>0&&fxx<0; t='max';
    else;               t='saddle'; endif
    fval = cx^3-3*cx+cy^3-3*cy;
    printf("  (%+d,%+d): D=%.1f → %s  f=%.1f\n", cx,cy,D,t,fval);
  endfor
endfor

% ─── Laplace's equation solver (finite differences) ──────────────
N=25; u=zeros(N);
u(end,:) = sin(linspace(0,pi,N));  % top BC: u=sin(πx)
for iter=1:3000
  u(2:end-1,2:end-1) = ...
    (u(1:end-2,2:end-1)+u(3:end,2:end-1)+ ...
     u(2:end-1,1:end-2)+u(2:end-1,3:end))/4;
endfor
% Exact: u(x,y) = sin(πx)·sinh(πy)/sinh(π)
exact_u = @(x,y) sin(pi*x).*sinh(pi*y)/sinh(pi);
x_pts = linspace(0,1,N); y_pts = x_pts;
[Xg,Yg] = meshgrid(x_pts, y_pts);
err_laplace = max(abs(u(:) - exact_u(Xg(:),Yg(:))));
printf("\nLaplace solver (25×25 grid) max error: %.3e\n", err_laplace);