Help me in implementing douglas-Peuker algorithm in matlab

Hi, i want to implement the "Douglas-Peuker" algorithm in matlab that is helpful to in my project in simplifying and determining the shapes of objects.

http://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm the above link has a simple explanation of the algorithm and a code demo, the thing i don't understand here is what is 'epsilon'. And there he used a recursion of the function, i don't understand how to implement recursion in matlab.

the code that i tried is

 function index= poly_approxd(binary)
%function to find approximate shape of objects, input here is a binary
%image. and the output is a list of points on boundary of objectswhich are
%results of Douglas-peuker algorithm.
bounds=bwboundaries(binary);
num_obj=size(bounds,1); %to find the number of objects, and to implement on all objects
for n= 1:num_obj
    list_pts=bounds{n};   %assigning points on boundary to a variable called list_pts.
    P1=[list_pts(1,2), list_pts(1,1)];   %starting point of boundary
    P2=[list_pts(end,2),list_pts(end,1)];  %last point on boundary
    num_pts=size(list_pts,1);  
    index=[];   % the list of points that that are to be output.
    for p=2:num_pts-1
        Pt=[list_pts(p,2), list_pts(p,1)];
        dmax=0;
        d=normal_dist(P1, P2, Pt);   %function to compute perpendicular distance of point to the line joining 1st(P1) and last(P2) points
        if d>dmax
            dmax=d;
        else
            index=[index;Pt];
            P1=Pt;
        end
    end
end
return

and the function normal_dist is

function d= normal_dist(P1, P2, P)
%this function computes the distance(d) of a point 'P'(x0,y0) from a line passing through
%two known points P1(x1,y1) and P2(x2,y2). it uses the formula
%d=asb((P2(1)-P1(1))*(P1(2)-P(2))-(P1(1)-P(1))*(P2(2)-P1(2)))/sqrt((P2(1)-P1(1))^2+(P2(2)-P1(2))^2)
%that is d=|(x2-x1)(y1-y0)-(x1-x0)(y2-y1)|/sqrt((x2-x1)^2+(y2-y1)^2).

d=abs((P2(1)-P1(1))*(P1(2)-P(2))-(P1(1)-P(1))*(P2(2)-P1(2)))/sqrt((P2(1)-P1(1))^2 +(P2(2)-P1(2))^2);
return

but the result is not exactly what i'm expecting, the output list is including only one row, that is row of staring point. the problem i found here is 'bwboundaries' is not not listing boundary points in sequence as the curve is, it is giving the list with a starting point on the first column of the boundary.

can any one please help me in solving the problem... thank you in advance.