<p>In this paper, a fractional-order corneal shape model involving the Caputo fractional derivative is investigated. The proposed model extends the classical corneal shape equation by incorporating fractional-order dynamics, thereby providing a flexible framework for describing nonlocal effects. By employing suitable fixed-point techniques, existence of solutions is established via Krasnoselskii’s fixed point theorem, while uniqueness is obtained through Banach’s contraction principle under appropriate assumptions on the nonlinear term. Furthermore, Ulam–Hyers stability of the proposed fractional model is analyzed and sufficient conditions ensuring stability are derived. To support the theoretical analysis, an L1 finite difference scheme based on the Caputo fractional derivative is developed for the numerical approximation of the model. The effectiveness of the proposed numerical method is demonstrated through representative computational results. In addition, three illustrative examples are presented to validate the existence, uniqueness, and Ulam–Hyers stability results. The obtained analytical and numerical findings confirm the consistency and applicability of the proposed fractional corneal shape model.</p>