How can I convert longitude and latitude in x and y coordinates on a Mollweide map?

I have to display the coordinates attached to the question on a Mollweide map. To do so I've tried implementing the Newton iterative method shown in the wikipedia page "Mollweide projection", like so:

R = 6.378145*(10^6);
axis([-2*R*sqrt(2) 2*R*sqrt(2) -R*sqrt(2) R*sqrt(2)]);
told = latitude(1);
t = told;
teta = nan(length(latitude), 1);
x = nan(length(longitude), 1);
y = nan(length(latitude), 1);
% Iterative cycle
for i = 1 : length(latitude)
while (abs(2*t + sin(2*t) - pi*sin(latitude(i))) >= 1e-6)
    tnew = told - (2*told+sin(2*told)-pi*sin(latitude(i)))/(2+2*cos(2*told));
    told = tnew;
    t = tnew;
end
    teta(i) = t;
    x(i) = (R*2*sqrt(2)/pi)*longitude(i)*cos(teta(i));
    y(i) = R*sqrt(2)*sin(teta(i));
end

x = x*pi/180;

Now this doesn't seem to work since the attached points, that can be seen by running this simple code:

figure(2);
hold on;
currIdx = 1;
for i = [find(abs(diff(longitude')) > 10), length(longitude)] 
    idx = currIdx:i;
    if length(idx) < 2
        style = '.r';
    else
        style = 'r';
    end
    plot(longitude(idx), latitude(idx), style);
    axis([-180 180 -90 90]);
    currIdx = i+1;
end

defenitely don't match with those obtained with the first piece of code. They can be visualized in the same way:

figure(3);
axis([-2*R*sqrt(2) 2*R*sqrt(2) -R*sqrt(2) R*sqrt(2)]);
hold on;
currIdx = 1;
for i = [find(abs(diff(x')) > 1.5e7-1e7), length(y)] 
    idx = currIdx:i;
    if length(idx) < 2
        style = '.r';
    else
        style = 'r';
    end
    plot(x(idx), y(idx), style);
    currIdx = i+1;
end

Why does this happen and how can I get the correct coordinate conversion?

[Christian]

I need to do the same thing. I hope you don't mind me asking before I get into it myself, but did you find any solution up to now?