You are not updating the angle in your simulation as the velocity vector is changing. Here is a modified version of your 'animate' function and the result.
t.append(t[i] + dt)
vx.append(vx[i] + dt * ax[i]) # Update the velocity
vy.append(vy[i] + dt * ay[i])
vel = np.sqrt(vx[-1] ** 2 + vy[-1] ** 2) # magnitude of velocity
Alright I've taken a closer look a the code. One first advice that I could give you is to try to understand how the code works. It is very easy to understand that the for-loop is a time loop, i.e. the system is advanced forward one time step at each iteration of the loop. The time marching equation uses a theta-scheme for the time discretization (theta=0.5 --...
Your function func(rho, gamma) takes two arguments while you pass only one (I think none). That is what your error tells you.
Recheck where you use func(rho, gamma) and make sure you pass two arguments.
I was able to do it this way though it may had been easier to plot it the way Federico has mentioned before: imagesc(A).
I declared an A matrix of zeros of mxn dimension and then loaded the ratings. I then applied HeatMap(A) which yields the image below.
A = importdata('u.data');
user_id = A(:, 1);
movie_id = A(:, 2);
rating = A(:, 3);
% Build matrix R
What happens as you add more points to your domain near $x=x_c$? Also, many (most, all(?)) definitions of the Heaviside step function give a value of the mean of the extremes at the transition point, but you initialize to zero and skip $x=x_c$ in your if-block. You should handle that case, too.
Also, why write your own when numpy supplies one?
Edited to ...
I answered a similar question in StackOverflow.
The main "trick" is to transform your grid before plotting. For example, using the following transformation
&t' = t\, ,\\
&x' = x e^t\, ,
and then use $(t', x')$ for your plot with contourf() or pcolormesh().
Following is a snippet showing the main idea.
import numpy as np