MATLAB using ode45 to animate wave function throws tolerance error

I am working on animating a wave function in MATLAB. So far, I have not reached to part of the animation yet since my ode45 throws an error. I am using de Korteweg-de Vries equations for this (which can be found here: https://en.wikipedia.org/wiki/Korteweg%E2%80%93de_Vries_equation).

Now my code is as follows:

function[t, N, u] = wave(u0, L, T, S, N)

%First we create a vector t which stores the time entries, a vector x which
%stores the location entries and an h which is the location step size for
%the Korteweg-De Vries function.
time = linspace(0,T,S);
x = linspace(-L,L,N);
h = x(2)-x(1);

U0=u0(x);

options=odeset('RelTol',1e-13,'AbsTol',1e-13);
[t,u] = ode45(@kdv,time,U0,options);

    function dudt = kdv(t,u)

    d2udx2 = ([u(N)-2*u(1)+u(2); diff(u,2); u(N-1)-2*u(N)+u(1)] ./ h^2);
    total = d2udx2 + 3.*u.^2;

    diff_total = diff(total);
    diff_total = [diff_total(end); diff_total(1:end)];

    dudt = -[diff_total(2)-diff_total(N-1); diff(diff_total,2); diff_total(N-1)+diff_total(2)] ./ (2*h);

    end
end

Now when I set f=@(x)(2/2).*(1./cosh(sqrt(2).*x./2).^2) and then call the function with [t,N,u]=wave(f, 20, 10, 200, 200); I get the following error:

Warning: Failure at t=6.520003e-02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.220446e-16) at time t. In ode45 (line 360) In wave (line 23)

Also the returned t and u have a size of 16x1 and 16x200 respectively, but I think they would expand to 200x1 and 200x200 if the warning did not occur.

Any hints on how to solve this would be much appreciated.

Solution

You did too much differences, the middle step of constructing diff_total is redundant.

d2udx2 and thus total are constructed correctly. You could shorten the first construction as
```
d2udx2 = diff([u(N);  u;  u(1) ], 2)./h^2`;
```
The next step constructing diff_total is redundant with the difference construction in dudt. The diff_total construction is lacking a division by h, for some unknown reason you have difference order 2 in the middle of dudt.
The outer x differentiation would be best done via convolution to get a central difference quotient as you have in the outer elements,
```
dudt = -conv([ total(N); total; total(1)], [1; 0; -1], 'valid')/(2*h)
```

All together this gives the derivatives procedure

function dudt = kdv(t,u,h)
    d2udx2 = diff([u(end);  u;  u(1) ], 2)./h^2;
    total = d2udx2 + 3*u.^2;
    dudt = -conv( [ total(end); total; total(1)], [1; 0; -1], 'valid')./(2*h);
end

which (with N=80) produces the solution

which as intended is a soliton wave traveling with speed c=2.