A = 2; B = 4;

\(eplus[{x1_, y1_}, {x2_, y2_}, p_]\  := \ 
      Module[
		{la, beta, inv1, inv2, xs, ys}, 
		
        If[\((x1 == 0)\) && \((y1 == 0)\), Return[{x2, y2}]]; 
		
        If[\((x2 == 0)\) && \((y2 == 0)\), Return[{x1, y1}]]; 
		
        If[\((x1 == x2)\) && \((Mod[y1 + y2, p] == 0)\), Return[{0, 0}]]; 

        		inv1\  = \ 
          Mod[\(\(ExtendedGCD[2\ y1, p]\)[\([2]\)]\)[\([1]\)], p]; 
        inv2\  = \ Mod[\(\(ExtendedGCD[x2 - x1, p]\)[\([2]\)]\)[\([1]\)], p]; 
        
		la\  = \ 
          Mod[If[x1 == x2, \((3\ x1^2\  + \ A)\)\ inv1, \((y2 - y1)\)\ inv2], 
            p]; 
		beta\  = \ Mod[y1 - \ la\ x1, p]; 
		
        xs\  = \ Mod[la^2\  - x1\  - x2, p]; 
		
        ys\  = \ Mod[\(-\((la\ xs\  + \ beta)\)\), p]; 
		{xs, ys}]\)

eplus[{4, 1}, {0, 3}, 5]

\(For[a = {0, 0}; b = {0, 2}; i = 0, i < 7, \(i++\), 
	
      Print[a = eplus[a, b, 5]]]\)